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ABSTRACT 


The  heat  checking  phenomenon  In  mechanical  seals  resulting  from 
the  passage  of  asperities  traversing  at  a  high  speed  over  the  surface 
of  a  seal  ring  is  solved  with  a  thermomechanical  model  of  a  seml- 
Inflnlte  medium  subjected  to  frictional  heating  over  a  moving  contact 
zone.  The  traverse  speed  is  considered  to  be  the  same  as  the  rubbing 
speed  for  the  current  problem.  The  two-dimensional  theory  depicts 
essentially  a  single  moving  line  contact  with  uniform  or  parabolic 
pressure  distribution  over  the  narrow  width.  For  the  actual  asperity, 
single  or  periodic,  the  contact  area  Is  best  described  by  a  “spot"  that 
may  be  approximated  with  a  circular  or  a  rectangular  geometry.  For 
such  a  configuration,  a  three-dimensional  theory  is  developed.  In  the 
present  analysis  the  general  solutions  are  obtained  for  arbitrary 
distribution  of  contact  loads.  Numerical  solutions  are  given  for  both 
uniform  and  non-uniform  pressure  distributions.  The  roles  of  material 
properties  and  operating  variables  are  expressed  In  terms  of  dimension¬ 
less  parameters.  The  double  Fourier  transform  and  the  Green's  function 
methods  are  used  throughout  the  analysis.  Numerical  schemes  for  the 
corresponding  analytic  solutions  of  integral  form  are  developed. 

For  a  constant  traversing  speed  of  15  m/s  (600  ips)  and  a  Coulomb 
coefficient  of  0.5,  the  two-dimensional  theory  for  a  single  moving 
asperity  concludes  that  the  crack  could  initiate  at  a  normal  average 
pressure  of  482  MPa  (70,000  psi }  for  the  case  of  a  parabolic  pressure 
distribution  over  the  narrow  width  of  1  mm  (0.04  in.).  The  three- 
dimensional  theory  for  a  single  moving  asperity  with  a  circular  contact 


area  of  radius  0.25  ram  (0.01  In.)  yields  a  critical  average  pressure  of 
196  MPa  (28,400  psl)  at  the  fracture  threshold. 

The  Investigations  show  that  the  distribution  of  the  pressure  over 
the  contact  zone  Is  more  of  a  deciding  factor  than  the  shape  of  the 
contact  area.  In  addition,  heat  checking  from  cyclic  asperity  excita¬ 
tion  Is  worse  than  that  from  a  single  asperity  excitation.  However, 
all  results  predict  a  be low- the- surface  crack  initiation.  Typically, 
for  an  asperity  with  the  half-width,  or  radius,  of  the  contact  area 
equal  to  0.5  mm  (0.02  In.),  the  critical  location  for  crack  Initiation 
is  0.05  mm  (0.002  In.)  below  the  surface.  The  result  in  a  sense 
justifies  the  use  of  linear  thermoelastic  theory,  since  at  this  depth 
the  material  considered  remains  In  an  elastic  state. 

The  cumulative  effect  Is  demonstrated  with  a  sequence  of  three 
asperities  spaced  apart  by  a  distance  of  \  «  12  (12  times  the  asper¬ 
ity  characteristic  dimension  a).  As  compared  with  196  MPa  (28,400  psl) 
for  a  single  moving  asperity,  the  critical  average  pressure  becomes  88 
MPa  (12,700  psi)  which  Is  much  lower.  Consequently,  the  possibility  of 
Initiating  cracks  Is  higher  when  periodic  moving  asperities  can  occur 
in  the  contact  zone  than  when  only  a  single  asperity  can  occur. 
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CHAPTER  1 


INTRODUCTION  AND  LITERATURE  SURVEY 


1.1  Purpose  of  the  Investigation 

The  Investigation  uses  several  mathematical  models  to  study  the 

possibility  of  fracture  initiation  by  dry  friction  and  the  configura¬ 
tion  of  fracture.  The  general  study  of  friction  cracking  has  a  special 
application  in  marine  seals.  The  purpose  of  this  study  is  to  extend 
our  understanding  of  the  problem  by  looking  Into  the  thermomechanical 
state  of  stress  in  a  mechanical  part  resulting  from  the  traversing  of  a 
single  asperity  or  periodic  asperities  over  its  surface.  The  under¬ 
standing  of  such  a  failure  process  could  be  used  to  improve  the  design 
of  devices  such  as  marine  seals  by  alleviating  the  problem  of  friction 
cracking. 

1.2  General  Background 

When  two  bodies  are  in  high  speed  sliding  contact  under  heavy 
loads,  the  localized  contact  area  is  much  smaller  than  the  design  slid¬ 
ing  surface.  The  resulting  pressure  in  the  area  can  ce  higher  than  the 
design  pressure  by  an  order  of  magnitude.  Localized  hot  spots  may 
occur  as  a  result  of  excessive  frictional  heating  near  the  contact  sur¬ 
faces.  Because  of  a  combination  of  thermal  heating  and  the  mechanical 
load,  the  material  may  crack  in  the  neighborhood  of  the  contact  zone. 
This  phenomenon  is  called  "heat  checking"  or  "thermocracking".  It  com¬ 
monly  occurs  in  mechanical  seals  and  brakes  as  shown  in  Figures  1.1  and 
1.2.  In  general  numerous  radial  cracks  can  be  observed  perpendicular 


Figure  1. 


.  Radial  hairline  cracks  on  the  metallic  ring 
after  running  against  a  carbon  ring  at  a  high 
peripheral  speed. 
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Figure  1.2.  Thermal  cracks  on  the  brake  shoe. 


to  the  sliding  direction  and  almost  periodically  separated  along  the 
circumference.  While  extensive  analysis  of  thermoelastic  effects  In 
brakes  has  been  studied  [1-5]1,  the  b^0.  of  which  Is  in  the  railroad 
Industry,  much  intensive  study  of  the  effects  is  needed  in  mechanical 
seals.  Of  concern  here  will  be  the  face  seals  located  along  the  pro¬ 
peller  shaft  in  submarines  to  prevent  leakage  of  sea  water  through  the 
shaft  tunnel.  Figure  1.3  illustrates  one  assembly  of  such  a  device,  in 
which  the  seal  consists  essentially  of  two  annular  rings,  one  of  which 
is  fixed  to  the  shaft  tunnel  housing,  another  Is  mounted  on  and  rotates 
with  the  shaft.  The  mating  surfaces  of  the  rings  therefore  are  pressed 
against  each  other  under  spring  pressure  and  rub  at  a  high  speed.  The 
design  nominal  pressure  between  the  seal  rings  depends  on  the  type  of 
vessel  and  the  location  of  the  seal.  For  Instance,  a  low  pressure  of 
about  100  kPa  (14.5  psl)  Is  sufficient  for  most  surface  vessels;  for 
submarines,  the  pressure  required  could  Increase  by  orders  of  magni¬ 
tude.  Had  the  pressure  been  evenly  distributed  according  to  design, 
the  service  life  of  the  face  seals  would  no  longer  be  a  serious  problem 
even  at  a  high  rubbing  speed.  However,  It  Is  well  known  that  the 
actual  contact  area  may  only  be  a  small  fraction  of  the  nominal  area  at 
the  design  Interface.  In  other  words,  a  low  nominal  design  pressure 
may  very  well  result  in  a  very  high  interfacial  pressure,  thus  a  very 
high  dry  frictional  force  in  the  actual  contact  area.  The  high  fric¬ 
tion  would  cause  locally  an  extremely  high  temperature,  called 

1  Numbers  in  brackets  denote  references. 


"flash  temperature"  by  Archard  [6].  The  local  contact  area  is  also 
called  the  “red  banding"  or  “hot  spot"  [7],  which  has  been  experimen¬ 
tally  demonstrated.  In  severe  cases  the  temperature  can  be  extremely 
high,  leading  to  cracking  of  the  surface  [8].  If  we  use  the  figure  of 
10"3  as  area  ratio  (contact  area/nominal  area) — Burton  [9]  considered 
10"4  as  a  possible  area  ratio— a  low  design  pressure  of  240  kPa  (35 
psl)  could  result  in  a  240  MPa  (35,000  psi)  local  pressure  In  the  con¬ 
tact  zone,  a  pressure  well  within  the  range  of  fracture  Initiation  with 
unfavorable  Coulomb  coefficients. 

The  cause  of  the  localization  of  the  contact  area,  which  is  In  the 
province  of  thermoelastic  instability  (TEI),  will  not  be  considered; 
nor  will  the  metallurgical  change  In  the  surface  layer.  It  will 
suffice  to  assume  that  the  localization  of  contact  area  is  due  to  some 
form  of  asperities,  which  may  be  fixed  to  any  of  the  mating  surfaces  or 
may  precess  with  respect  to  both.  The  present  work  treats  the  thermo- 
mechanical  stress  state  in  the  part  that  mates  with  the  one  containing 
the  asperity  and  the  possibility  of  heat  checking  therein.  The  speci¬ 
fic  investigation  will  consider  the  asperity  or  asperities  being  on  the 
carbon-graphite,  which  moves  relatively  on  the  face  seal  of  Stellite 
III.  The  failure  due  to  heat  checking  is  on  the  face  seal  which  prior 
to  the  onset  of  failure  is  assumed  to  have  no  other  defect  on  its  con¬ 
tacting  surface. 

1.3  Related  Investigation  in  Progress 

The  failure  records  have  shown  that  the  performance  life  of 
devices  with  pressurized  moving  contact  surfaces  becomes  drastically 
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reduced  whenever  the  parameters  such  as  the  load,  sliding  speed,  or  the 
temperature  field  attain  certain  range  of  values.  There  has  been 
Increased  emphasis  on  finding  a  solution  of  failure  control,  experimen¬ 
tally  and  analytically,  of  this  physical  problem  In  recent  years.  A 
general  survey  of  the  problem  of  cracking  through  the  development  of  a 
frictional  hot  spot  was  discussed  by  Burton  [7].  A  series  of  experi¬ 
ments  carried  out  by  Sibley  and  Allen  [10]  showed  photographic  evidence 
of  systematically  moving  hot  patches  In  the  contact  zone.  Based  on  the 
results  In  the  recent  literature  [11-15],  It  can  be  assumed  that  TEI 
precedes  thermocracking  and  that  contact  Is  localized  at  several  hot 
patches.  Considerable  evidence  was  also  shown  by  Kennedy  [16]  to  sup¬ 
port  this  crucial  assumption.  In  his  experiment,  using  the  carbon  ring 
against  a  metallic  mating  ring  made  from  440  C  stainless  steel,  beryl¬ 
lium  copper  or  52100  bearing  steel  under  both  dry  and  liquid  lubricated 
conditions,  the  existence  of  distinct  spot  asperities  on  the  metallic 
ring  was  observed.  It  was  found  that  the  spots  at  least  In  dry  opera¬ 
tion  tend  to  remain  stationary  with  respect  to  the  metallic  mating  ring 
of  the  seal,  whether  that  ring  be  stationary  or  rotating.  However, 
other  Investigations  have  shown  hot  patches  moving  relative  to  the  mat¬ 
ing  ring  and  stationary  on  the  primary  ring  [12,13].  In  addition. 
Burton  [14]  reported  that  for  an  aluminum  ring  sliding  on  a  glass  disk, 
the  hot  spot  precessed  at  a  much  lower  speed  than  the  rubbing  speed. 

The  uncertainty  of  this  observed  discrepancy  on  the  speed  of  the  moving 
asperities  remains,  but  there  Is  no  doubt  about  the  existence  of  the 
moving  asperities  due  to  TEI  on  mechanical  face  seals  through  the 
effort  of  all  the  researchers  in  this  field.  Therefore,  several 
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analytical  studies  of  the  failure  due  to  the  existence  of  the  moving 
asperities,  which  are  characterized  by  high  normal  and  frictional  loads 
and  the  resulting  thermal  heating  on  the  contact  surface  of  the  sliding 
mechanical  parts,  have  been  In  progress. 

When  the  thermal  effect  Is  neglected,  the  study  Is  isothermal. 
Various  analyses  of  stresses  and  displacements  Induced  In  an  elastic 
half  plane  by  a  concentrated  line  load  moving  at  a  constant  speed  along 
Its  surface  have  been  developed  [17*19].  Eason  [21]  dealt  with  several 
types  of  moving  surface  forces  and  procured  the  solution  of  the  mechan¬ 
ical  stress  and  displacement  field  In  a  semi-infinite  solid.  Surface 
displacement  and  temperature  field  of  a  convective  elastic  half  space 
under  an  arbitrarily  distributed  fast-moving  line  heat  source  were 
obtained,  using  Integral  transform  methods,  by  Ling  and  Mow  [20].  The 
problem  of  thermal  stresses  [22]  and  temperature  distribution  [23]  was 
examined.  A  finite  element  analysis  was  recently  developed  by  Kennedy 
[24]  for  studying  surface  temperatures  resulting  from  frictional  heat¬ 
ing  in  sliding  systems.  He  also  applied  finite  element  techniques  to 
study  the  stresses  In  the  mechanical  face  seal  [25]  and  showed  that  the 
dominant  stresses  in  the  seal  components  are  thermal  stresses. 

Marscher  [26,27]  analyzed  the  thermomechanical  stress  state  in  the  part 
that  contains  the  asperity.  A  two-dimensional  repeated  loading  phenom¬ 
enon  was  first  studied  for  stability  of  contacting  surfaces  by  Burton, 
et  al .  [28,29].  The  surface  stress  component  (parallel  to  the  surface) 
resulting  from  a  periodic  row  of  moving  patches,  with  width  Zi  each, 
and  a  spacing  of  2m  was  investigated  by  Tseng  and  Burton  [30],  who 
restricted  the  analysis  to  the  two-dimensional  plane  stress  case.  They 


concluded  from  the  computed  results  that  the  tensile  stress  would 
appear  Instantaneously  with  each  passage  of  the  heat  source. 

1.4  Mathematical  Model  and  Analytical  Approach 

According  to  the  results  of  the  experiment  by  Kennedy,  et  al. 

[16],  the  width  of  the  moving  asperity  was  about  0.1  to  1  rid  (0.004  to 
0.04  in.);  while  the  diameter  of  a  typical  marine  shaft  seal  could  be 
50  an  (19.7  in.)  or  more.  Because  of  this  size  difference  between  the 
contact  area  and  the  seal,  as  shown  in  Figure  1.4,  the  mathematical 
model  Is  represented  by  a  half-space  subjected  to  fast-moving  loads  of 
a  constant  velocity  (Y).  The  loads  come  from  a  combination  of  thermal 
and  traction  contributions.  Fast  moving  refers  to  the  case  wherein  V 
»  </a,  where  <  and  a  are  defined  as  the  thermal  dlffuslvlty  and  the 
half  width  of  the  asperity.  The  Investigation  starts  with  the  study  of 
a  model  with  a  moving  line  asperity.  For  such  an  asperity  a  two- 
dimensional  analysis  Is  adequate  to  develop  the  thermomechanical  solu¬ 
tions.  In  order  to  better  model  the  contact  area  and  periodic  excita¬ 
tion  of  asperities  in  practice,  the  research  adopts  eventually  a  more 
realistic  three-dimensional  analysis.  Since  It  has  been  observed  that 
the  lubricated  film  Is  continually  being  broken  by  small  patches  of 
solid-to-solid  contact  during  seal  operation  [16],  we  may  as  well 
assume  that  the  working  seals  are  In  dry  contact  representing  the  worst 
operating  condition.  In  addition,  for  the  worst  case,  the  mating  seal 
that  contains  the  asperity  is  regarded  as  an  insulator.  The  analyses 
are  based  on  the  general  theory  of  a  continuum.  The  relative  speed  of 
contact  surfaces  is  assumed  to  be  large  enough  to  result  in  a  high 
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Peclet  number  (Ya/<).  But  it  Is  much  smaller  than  the  Rayleigh  wave, 
which  for  a  steel  based  seal  material  is  approximately  2800  m/s  (11  x 
104  ips),  and  justifies  the  use  of  the  guasl-static  theory.  Further¬ 
more,  since  such  loading  also  results  in  a  low  dilatation  rate,  the 
uncoupled  thermal  stress  theory  is  then  applicable.  From  observations 
of  failed  specimens  near  the  surface,  it  is  reasonable  to  postulate 
that  the  plastic  deformation  and  rupture  at  the  surface  are  at  the 
granular  or  even  the  subgranular  level.  The  base  solid  material  sub¬ 
jected  to  the  asperity  friction  is  essentially  elastic.  Consequently, 
the  linear  theory  of  thermoelasticity  holds.  The  principal  advantage 
of  using  the  linear  theory  is  the  application  of  the  superposition 
principle,  enabling  a  delineation  in  the  stress  field  to  one  from  the 
mechanical  load  of  the  moving  pressure  and  friction  and  another  from 
the  moving  heat  input,  a  result  of  the  irrecoverable  frictional  work. 
The  combined  effects  will  then  determine  the  possibility  of  fracture 
initiation  using  the  tensile  fracture  criterion.  In  two-dimensional 
analysis,  the  Green's  function  approach  is  used;  while  in  three- 
dimensional  analysis,  the  analytical  solutions  of  the  mechanical  and 
the  thermal  stress  fields  are  solved  by  means  of  double  Fourier  trans¬ 
form.  The  temperature  field  is  obtained  by  both  the  Integral  transform 
and  Green's  function  methods.  Numerical  results  of  the  corresponding 
Integral  solutions  are  represented  in  graphs.  Schematic  geometries  and 
detailed  analyses  are  given  in  the  following  chapters. 

Two-dimensional  and  three-dimensional  models  are  examined  in 
Chapters  2  and  3,  respectively.  The  former  demonstrates  a  simpler 
means  of  formulation;  while  the  latter  gives  a  better  description  of 
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the  physical  problem.  Both  models  deal  with  a  single  moving  asperity. 
In  Chapter  4,  a  three-dimensional  model  of  periodic  moving  asperities 
is  investigated  further.  The  formulation  leads  to  a  general  solution 
which  can  be  used  to  justify  the  results  from  Chapters  1  and  2. 
Detailed  discussions  and  conclusions  are  given  in  Chapter  5  and  the 
related  mathematical  background  is  given  in  the  appendices. 


1.5  The  Uncoupled,  Quasi-Static  Thermoelastic  Theory 

The  analysis  is  restricted  to  the  dynamic  and  coupled  theory  of  an 
isotropic  linear  thermoelastic  solid,  which  can  be  found  in  most 
thermoelasticity  textbooks  [31-36].  In  vector  form,  the  displacement 
formulation  in  the  absence  of  body  forces  and  the  Duhamel  heat  equation 
[37]  with  no  internal  heat  generation  are 

u^u,  +  (x  +  u)  grad  div  u  -  (3x  +  2U)  a  grad  T  *  pu  ,  (1.1) 


kv2T 


pcT 


1  + 


(3X  +  2y)  a  T  e 


pc 


o  "kk 
aT 


(1.2) 


where  T0  is  the  temperature  of  the  solid  in  the  unstressed  state, 

.  • 

eick  Is  the  dilatation  rate  and  T  is  the  temperature  rate.  Equations 
(1.1)  and  (1.2)  are  simplified  by  deleting  the  inertia  term,  0U  for  the 
quasi-static  case  and  of  the  mechanical  coupling  term. 


(3x  +  2ji)  a  Tq 


pc 


aT 


for  the  uncoupled  theory.  Thus  the  governing  equations  for  the 
uncoupled,  quasi-static  thermoelastic  theory  are 


uV  u  +  (x  +  y)  grad  div  u  *  (3x  +  2p)  a  grad  T 


(1.3) 


72T  *  -  T  .  (1.4) 

K 

A  short  discussion  follows  to  justify  the  use  of  the  quasi-static 
uncoupled  theory.  The  mechanical  coupling  term  is  considered  first. 

It  can  be  shown  that  for  steel  it  Is  approximately 


1  +  10*5  T 

0  «T  J 

where  T  is  absolute  temperature  measured  in  Rankines.  For  aluminum,  it 
becomes 


kV2T  *  pcT 


kV2T  *  pcT 


1  +  1.8  x  10 


If  SI  units  are  used,  l.e.,  the  temperature  Is  measured  in  degrees 
Kelvin,  the  coefficients,  10“ 5  for  steel  and  1.8  x  10" 5  for  alu¬ 
minum  should  be  replaced  by  1.8  x  10"5  and  3.24  x  10"5,  respective¬ 
ly.  It  is  argued  that,  for  most  cases  of  a  high  thermal  load,  ejck/ 

<xf  should  be  of  the  order,  0(1),  so  that  the  coupling  term  is  negli¬ 
gible  except  for  conditions  in  which  the  temperature  distributions  have 
sharp  variations  or  discontinuities  in  their  time  histories,  which 
often  occurs  during  the  propagation  of  thermoelastic  waves  in  the 
aftermath  of  thermal  shocks.  The  similar  argument  is  also  made  by 
Boley  and  Weiner  [38]  who  rewrote  Equation  (1.2)  as 
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kv  T  »  pcT 


L  aT  J 

where  $  Is  defined  by 

(3x  +  2W)2  a2TQ 
^  (X  +  2y)  pC 

which  has  the  values  0.029  and  0.014  for  aluminum  and  steel  at  T0  * 
660*R  (corresponding  to  367  K),  respectively.  Further,  they  showed 
that  for  both  cases  the  coupling  Is  small  If 

• 

ekk 

«  20  . 

3  aT 

They  therefore  concluded  that  the  possibility  of  omitting  the  coupling 
terms  depends  not  only  on  the  fact  that  the  Inequality  6  «  1  must  hold 
(as  It  does  for  most  metals),  but  also  on  the  fact  that  strain  rates 
must  be  of  the  same  order  of  magnitude  as  temperature  rates  times  the 
coefficient  of  thermal  expansion. 

The  dynamic  effect  can  be  a  result  from  either  a  dynamic  loading 
state  or  a  non-steady  thermal  state  in  which  the  time  rate  of  tempera¬ 
ture  change  could  keep  up  with  the  stress  waves  in  the  material.  It 
was  stated  by  Ouhamel  [39]  that  the  inertia  term  can  be  disregarded  if 
the  time  rate  of  change  of  temperature  is  slow  enough.  Parkus  [40] 
showed  that  the  significant  effect  from  inertia  term  can  arise  only 
when  there  Is  an  Instantaneous  change  In  the  surface  temperature  or  In 
the  temperature  of  the  surrounding  medium.  In  fact,  the  dynamic  effect 
Is  greatly  reduced  if  the  temperature  change  occurs  in  a  very  short, 
but  finite,  interval  of  time.  This  was  confirmed  by  Danilovskaya 


[41-42],  who  studied  the  dynamic  effect  due  to  a  thermal  shock  on  the 
surface  of  a  half-space  in  detail  and  demonstrated  that  the  maximum 
dynamic  stress  is  reduced  to  86  percent  even  for  the  extremely  short 
duration  of  heating  of  10"12  seconds.  In  general,  under  usual  condi¬ 
tions  of  heat  exchange,  the  rate  of  temperature  change  is  small  in  com¬ 
parison  with  the  speed  of  sound  in  the  material.  Thus,  at  any  instant, 
the  thermal  stress  state  can  be  determined  by  the  instantaneous  values 
of  the  temperature  field.  Consequently,  thers  is  no  need  to  consider 
the  inertia  force  corresponding  to  the  motion  of  the  particles  during 
the  varying  thermal  expansion,  i.e.,  the  Inertia  term  p£  in  Equation 
(1.1)  can  be  omitted. 

1.6  The  Brittle  Fracture  Criteria 

Metals  can  be  generally  divided  Into  two  categories:  (1)  ductile 
metals.  In  which  marked  plastic  deformation  commences  at  a  fairly 
definite  stress  and  which  exhibits  considerable  ultimate  elongation 
prior  to  fracture;  and  (2)  brittle  metals,  for  which  plastic  deforma¬ 
tion  is  not  clearly  evident  at  fracture.  Since  the  heat  checking 
phenomenon  as  evidenced  in  Stellite  III  and  in  ceramic  materials  [43] 
is  associated  with  very  little  plastic  flow,  the  brittle  fracture  cri¬ 
terion  is  ascertained.  Based  on  the  hypothesis  enunciated  by  Griffith 
[44]  in  1921,  the  weakening  of  the  material  is  attributed  to  the  growth 
of  minute  cracks  resulting  from  an  induced  stress  field.  He  [45] 
derived  the  criterion  for  failure  under  a  biaxial  principal  stress 
field  <j\  and  a 3  at  an  arbitrary  orientation  8  with  respect  to  the  axis 
oj  as  seen  in  Figure  1.5,  and  studied  in  detail  the  variation  of  the 
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peripheral  stress  ?  on  the  surface  of  a  flat  elliptical  crack  (5  * 
5o).  He  postulated  that  fracture  would  result  when  the  maximum  periph¬ 
eral  stress  reaches  a  characteristic  value  of  the  material.  His  theory 
modified  by  McClIntock  and  Walsh  [72],  who  took  into  account  the  fact 
that  cracks  may  be  expected  to  close  under  sufficiently  high  compres¬ 
sive  stresses,  is  stated  as  follows.  As  the  normal  stress  <jn  acting 
across  the  face  of  a  pre-existing  crack  causes  the  crack  to  close, 

I.e.,  when 


on  3  C(o^  +  Og)  +  (<j^  -  03)  cos  2s]  <  0 


or 


ojCUf  +  1)*^  +  iif]  + 


o3[(4  +  u1/2 


Ilf]  <  0  > 


then  a  compressive  fracture  [72]  will  occur  when 


<4 +  D1/2  +  vf 
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where  ac  is  the  uniaxial  compressive  strength  of  the  material,  is 
the  surface  coefficient  of  friction  and  9,  which  is  equal  to  0.5 
tan"1  (i/uf)  and  measured  from  the  minor  principal  stress  <j3  axis,  is 
the  orientation  of  the  crack  at  which  fracture  initiation  will  com- 
mer  -  ,  If  0,  9  ♦  tt/4;  if  +  «,  9  +  0.  Thus  the  direction 

of  fracture  is  always  inclined  at  an  acute  angle  to  the  direction  of 
the  axis  of  the  largest  principal  compressive  stress. 

When  an  causes  the  crack  to  open,  i.e.,  as  <?n  >  0,  tension  crack¬ 
ing  occurs.  For  the  case  that 


fracture  initiation  occurs  when  oi  *  au  where  <yu  is  the  uniaxial 
tensile  strength  of  the  material.  For  this  condition  of  tensile  fail¬ 
ure,  the  crack  begins  to  extend  In  its  own  plane  in  a  direction  perpen¬ 
dicular  to  that  of  the  major  principal  stress  <jx.  While  for  the  case 
that 

°3  <  ‘3  °u  * 

failure  takes  place  if 

(oi  '  a3)Z  +  8  °uUl  +  c3)  *  0  » 

and  the  crack  runs  towards  the  o3  axis  as  in  compressive  fracture.  The 
corresponding  fracture  locus,  as  shown  in  Figure  1.6,  is  well  illus-- 
trated  by  McClintock  [46].  He  also  extended  the  criterion  to  the 
three-dimensional  case  (see  Figure  1.7). 

Concerning  the  crack  propagation  direction  under  uniaxial  compres¬ 
sion  of  plates  of  brittle  materials,  both  analytical  and  experimental 
results  from  Nemat-Nasser,  et  al .  [73]  provide  good  correlation  with 
the  Griffith's  theory.  It  was  shown  that,  for  the  wide  range  of  pre¬ 
existing  crack  orientation,  the  out-of-plane  crack  extension  initiates 
at  an  angle  close  to  72*  from  the  direction  of  the  pre-existing  crack, 
and  then  curves  into  a  position  almost  parallel  to  the  maximum  axial 
compressive  stress  as  shown  in  Figure  1.8. 

Within  the  framework  of  the  linear  elastic  fracture  mechanics 
(LEFM),  three  modes  of  cracking  as  in  Figure  1.9  are  often  referred. 


Figure  1.5.  Brittle  fracture  under  biaxial  stress. 


Figure  1.7.  Brittle  fracture  solid  representing  the  states  of 
triaxial  stress  for  brittle  fracture. 


Of  concern  here  Is  the  Mode  I  crack  since  It  Is  often  observed  on  the 
failed  seal  [16].  Consider  a  through- the- thickness  microcrack  of 
length  2a  in  an  Infinite  plate  which  Is  subjected  to  the  uniaxial 
ultimate  tensile  stress  as  shown  in  Figure  1.10.  The  maximum  stress 
occurs  near  the  crack  tip  at  e  a  0  and  has  the  value  of  cr9d  * 
au/a/2r.  Under  this  condition,  the  Instability  of  the  microcrack 
results  In  the  fracture  initiation  predicted  by  Griffith's  theory. 


Figure  1.9.  The  three  modes  of  cracking. 
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Figure  1.10.  Crack  in  an  infinite  plate 


CHAPTER  2 


A  TWO-DIMENSIONAL  MODEL  OF  A  SINGLE  MOVING  ASPERITY 

2.1  Mathematical  Model2 

In  the  present  model,  it  is  reiterated  that  the  line  asperity  is 
on  the  insulated  mating  part  of  the  cracked  seal  and  the  speed  of  the 
effective  load  is  subsonic.  The  geometry  of  the  contact  area,  as  shown 
in  Figure  2.1  is  such  a  long  narrow  strip  that  the  width  of  the  area  in 
the  direction  of  motion  is  of  small  order  to  the  radius  of  the  seal 
surface.  Accordingly,  the  latter  may  be  represented  by  a  planar  half¬ 
space.  The  problem  may  then  be  described  by  a  plane  thermoelastlc 
model  as  shown  in  Figure  2.2.  The  contact  load  Is  represented  therein 
by  a  symmetrically  distributed  pressure  over  the  contact  zone.  Coulomb 
friction  is  present,  resisting  the  relative  motion  of  the  contacting 
surfaces,  that  is,  the  friction  force  is  in  the  traverse  direction  of 
the  load.  Two  sets  of  coordinates  are  employed:  (x^  -  x3)  are  fixed  to 
the  seal;  (xL  -  x3)  are  fixed  to  the  moving  load.  The  boundary  values 
of  this  problem  are  therefore  as  follows:  on  the  surface  boundary  x3  * 
0,  there  is  a  combined  normal  pressure  p(xL)  and  tangential  friction 
ufp(xl)  distributed  over  the  contact  zone  (-a  £  xj  a)  and  the 
half-plane  is  traction-free  elsewhere.  For  the  regions  at  infinity, 
regularity  conditions  hold. 

The  thermal  load  on  the  half-plane  comes  from  the  frictional  heat¬ 
ing  in  the  contact  zone.  The  mating  surface  is  by  postulation  adia¬ 
batic.  Therefore,  all  the  dissipative  energy  is  accounted  for  by  the 


2The  present  chapter  is  a  revision  of  [71], 


heat  Input  to  the  seal.  The  heat  loss  across  the  free  boundary  is  also 
neglected. 


The  material  properties  affecting  the  thermomechanical  stress 
field  will  be  the  elastic  coefficients  (E,v),  the  density  of  the  seal 
(p),  the  coefficient  of  friction  (uf)»  the  coefficient  of  thermal 
expansion  («),  the  thermal  conductivity  (It)  and  the  thermal  dlffuslvlty 
(k).  The  parameters  characterizing  the  loading  are  the  pressure  dis¬ 
tribution  p(x i )  over  the  contact  zone  and  the  velocity  of  the  moving 
load  (V)  for  the  case  when  the  traversing  speed  is  the  same  as  the 
rubbing  speed. 

The  stress  field  from  the  mechanical  loads  of  normal  pressure  and 
friction  will  at  first  be  computed.  The  thermal  load  from  the  friction 
force  will  then  be  considered  to  compute  the  temperature  rise  and  thus 
the  resulting  thermal  stress  field.  The  two  sets  of  stresses  should 
provide  an  insight  to  the  separate  effects  on  stresses  from  the  mechan¬ 
ical  and  the  thermal  loads.  The  combined  effect  from  superposing  the 
two  shall  provide  the  information  on  impending  failure  by  fracture. 

2.2  Mechanical  Stress  Field 

The  solutions  of  the  stress  fields,  corresponding  to  the  mechan¬ 
ical  loads  as  shown  in  Figure  2.2,  may  be  obtained  by  the  method  of 
impulse  response.  For  that  method  the  stress  field  {Sa^j}  from  3 
moving  normal  force  of  one  unit  (Figure  2.3)  is  the  impulse  response 
for  the  normal  pressure  while  from  3  moving  tangential  force  of 

one  unit  (Figure  2.4)  is  the  impulse  response  for  the  friction.  For 
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and  Cd  and  Cs  are  respectively  the  dilatational  and  the  shear  wave 
speeds. 

The  actual  stress  field  due  to  the  normal  pressure  Is  obtained  by 
superposition  as  the  Faltung  product  of  the  impulse  response  with  the 
pressure  distribution  function  p(xx),  i.e., 

°1j  *  5oij  *  p(xl)  *  <2.3] 

Similarly,  the  stress  field  due  to  the  friction  force  is  the  Faltung 
product  of  its  impulse  response  with  the  friction  distribution  func¬ 
tion: 
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We  shall  use  dimensionless  variables  such  that 
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where  A 


When  the  pressure  distribution  function  p(x)  is  given,  the  convo¬ 
lution  integral  can  be  evaluated.  The  simplest  case  Is  when  the  pres¬ 
sure  distribution  is  uniform,  “p  »  1,  for  which  Equations  (2.5)  and 
(2.6)  become 
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The  second  case  for  which  a  parabolic  pressure  distribution,  p  *  3/2(1 
-  X2),  is  considered  yields 
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2.3  Temperature  Field 

The  heating  on  the  boundary  Is  assumed  to  be  caused  principally  by 
the  friction.  The  analysis  is  therefore  one  with  a  half-space  sub¬ 
jected  to  a  fast-moving  heat  source.  When  the  mating  surface  and  its 
asperity  are  assumed  to  be  adiabatic,  the  rate  of  energy  dissipation 
from  the  friction  force  is  converted  entirely  into  the  heat  input  q 
such  that 


q  *  Vyfp(Xj)  . 

Using  the  solution  developed  by  Ling  and  Mow  [20],  we  may  obtain  the 
temperature  field  as 

♦(5.5)  *  0  ,  -«  <  e  <  -1  , 


UP. 


v/.  I 


I  V2  ) 


5  expj'  Ttf=Tr j dx  •  -l  <  5  < 1  • 


♦(C.c) 


UP. 


wf-l 


Q  ( X) 


p-,2  ) 


expf*7^n-j 


dx  ,  1  <  5  < 


(2.9) 


where 


♦U.c) 


Tk 


V 


is  the  dimensionless  temperature. 


Q(x) 


•2-  is  the  dimensionless  heat  input, 
qo 


P. 


Va 


is  the  Peclet  number. 


T  Is  the  temperature  above  the  stress-free  ambient  temperature  and  q0 
is  the  characteristic  heat  input  Uf%h  Q(x)  is  therefore  identified 
with  p(A). 

Surface  temperatures  for  the  case  of  uniform  and  parabolic  pres¬ 
sure  distributions  can  be  readily  obtained  by  letting  c  a  0.  The 
results  can  be  expressed  as 

♦s(5)  *  0  »  —  £  5  £  *1  » 
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respectively. 


2.4  Thermal  Stress  Field 

The  Impulse  response  of  a  stress  field  corresponding  to  the  heat 
Input  of  a  unit  pressure  Is  {5a},  which  Is  conveniently  separated  Into 
the  particular  and  the  complementary  solutions  {5cP}  and  {5oc}, 
respectively.  Their  expressions,  following  Mow  and  Cheng  [22] 3 ,  are 
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3The  equations  listed  are  a  corrected  version  of  those  given  In  the 
paper. 
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The  thermal  stresses  are  thus  computed  from  the  Faltung  product: 


°1j 


6 <Jij-  *  Q(  X) 


where  <Soij  Is  the  combination  of  the  particular  and  the  complementary 
solutions.  The  results  can  be  expressed  in  terms  of  convolution  inte¬ 
grals,  which  are 

ffij  a  0  ’  -  <  5  <  *1 


o-m 


/  p( X)  -  X,?)  dX  , 


-1  <  f  <  1 


.  *  j  pU)  So.-jU  -  X,?)  dx  ,  1  <  c  <• 


(2.14) 


in  which 
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Note  that  the  expressions  for  and  $3^.  correspond  to  those  in 
Equations  (2.12)  and  (2.13),  respectively,  with  the  deletion  of  the 
Heaviside  function. 

With  the  notation 
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Equation  (2.14)  for  the  uniform  pressure  loading  case,  on  the  sub¬ 


stitution  of  Equations  (2.15)  and  (2.16),  gives 
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for  1  <  5  <  «,  where 


Similarly,  the  stresses  for  the  case  of  the  parabolic  pressure 


f  -  X,c)  dx 


distribution  can  be  obtained  as  follows: 
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2.5  Numerical  Results 

The  numerical  evaluation  of  the  mechanical  stress  fields  in  Equa¬ 
tions  (2.7)  and  (2.8)  can  be  readily  obtained  by  direct  substitution  of 
the  data  for  loading  and  material  properties.  Figure  2.5  shows  * 
al  +  Op  for  different  depths  The  traversing  speed  V  *  15  m/s  (600 
ips).  The  solid  lines  represent  the  uniform  pressure  distribution  and 
the  dotted  lines  represent  the  parabolic  distribution.  Because  of  the 
load  discontinuity,  the  stress  component  ogg  has  two  singularities  on 
the  surface,  for  the  case  of  uniform  pressure,  at  {■  *  ±1.  In  the  real 
problem,  the  pressure  distribution  is  closer  to  the  case  of  a  parabolic 
function,  for  which  there  is  no  singularity. 

The  temperature  field  for  both  cases  is  shown  in  Figure  2.6.  Near 
the  surface,  there  exists  a  high- temperature  thermal  layer.  The  tem¬ 
perature  decays  exponentially  with  respect  to  the  depth  as  we  can 
observe  from  Equations  (2.9)  and  (2.10).  As  a  result,  the  thermal 
stresses  under  the  moving  load  may  be  compressive  near  the  surface  and 
become  tensile  at  a  certain  depth  where  the  maximum  temperature  gra¬ 
dient  occurs.  Indeed,  this  Is  the  case  and  illustrated  by  the  numeri¬ 
cal  evaluation  of  the  thermal  stress  fields  in  Equations  (2.17)  and 
(2.18),  and  shown  in  Figure  2.7. 

The  maximum  principal  stress  resulting  from  the  combined  stress 
field  is  shown  in  Figure  2.3  for  depths  5  a  1.0  and  0.1.  The 
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compressive  state  of  the  thermal  stress  In  the  surface  layer  (Figure 
2.7)  is  at  a  much  higher  level  than  that  from  the  mechanical  loading 
(Figure  2.5).  However,  the  combined  stress  field  does  become  tensile 
at  a  small  depth  and  reaches  the  maximun  principal  stress  at  approxi¬ 
mately  a  depth  of  one- twentieth  of  the  asperity  size.  According  to 
Figure  2.8,  a  tensile  crack  will  not  o^cur  because  (a.)  a  0.82 

is  less  than  one.  However,  the  temperature  and  the  combined  stress 
field  strongly  depend  on  the  asperity  size  and  the  other  parameters 
mentioned  before.  The  details  on  the  study  of  parameter  effects  will 
be  given  in  Chapter  5.  For  the  purpose  of  illustration,  the  maximum 
principal  stress  field  for  a  typical  asperity  size  of  1  mm  (0.04  in.) 
is  also  obtained.  These  results  are  compared  with  the  results  obtained 
with  an  asperity  size  of  0.5  nm  (0.02  In.)  in  Figure  2.9.  At  a  pres¬ 
sure  of  482  MPa  (70,000  psi ) ,  the  maximum  principal  stress  reaches  the 
ultimate  tensile  strength.  The  distributive  friction  force  of  241  MPa 
(35,000  psi),  using  a  coefficient  of  friction  =  0.5,  is  well  below 
the  ultimate  shearing  strength.  Therefore,  the  crack  could  initiate  at 
a  depth  of  0.05  mm  (0.002  in.)  from  the  surface  for  a  typical  asperity 
size  of  1  mm  (0.04  in. ) . 


Figure  2.9.  Principal  stress  fields  for  different  half  widths  (t 
0.25  mm  (0.01  in.),  l0  *  0.5  mm  (0.02  in.)). 


CHAPTER  3 


A  THREE-DIMENSIONAL  MODEL  OF  A  SINGLE  MOVING  ASPERITY 

3.1  Mathematical  Model 

In  order  to  model  a  hot  spot  with  a  contact  area  that  represents 
better  the  actual  case  and  to  investigate  the  three-dimensional  effect, 
the  three-dimensional  problem  of  a  single  asperity  moving  with  uniform 
velocity  over  the  surface  of  a  half  space  is  presented  in  this  chap¬ 
ter.  The  results  can  be  used  to  compare  with  the  two-dimensional 
results  obtained  previously.  This  half  space  representation  of  the 
seal  is  based  on  the  same  reasoning  concerning  the  relative  size 
between  the  seal  material  and  the  moving  asperity  as  used  In  the  two- 
dimensional  analysis.  The  asperity  contacts  with  the  material  on  the 
otherwise  traction-free  surface  in  a  small  rectangular  or  circular 
region.  The  coordinates  {xi }  are  fixed  to  the  moving  asperity  such 
that  xx  points  toward  the  trailing  direction  of  the  motion,  x2  is  per¬ 
pendicular  to  the  traversing  direction,  and  x3  Is  a  depth  measure 
pointing  from  the  surface  into  the  material.  Figure  3.1.  For  the 
rectangular  contact  area,  the  aspect  ratio,  t  *  b/a,  is  a  significant 
parameter  for  the  three-dimensional  effect.  For  large  aspect  ratio, 
the  asperity  excitation  is  effectively  a  moving  line  load  that  was 
presented  as  a  two-dimensional  plane  strain  solution  in  Chapter  2. 

In  the  three-dimensional  formulation  the  potential  theory  approach 
in  Chapter  2  is  no  longer  applicable.  Hence,  the  Fourier  transform 
method  (Appendix  I)  is  used  throughout  in  analyzing  the  mechanical  and 
the  thermal  stress  fields  as  well  as  the  temperature  field.  The  latter 


Is  also  analytically  checked  with  the  Green's  function  approach. 

Because  of  the  complexity,  the  general  solutions  given  in  Sections  2, 

3,  and  4  for  the  mechanical,  the  temperature,  and  the  thermal  stress 
fields  respectively  are  left  In  the  transformed  space.  The  numerical 
solutions  In  Section  5  apply  specifically  to  three  cases:  a  moving 
rectangular  contact  area  with  uniform  pressure  (Case  1),  a  moving  disk 
with  uniform  (Case  2),  or  nonuniform  (Case  3)  pressure.  In  the  Inverse 
Fourier  transform  the  multiple  Integrals  for  the  mechanical  stress 
field  are  reduced  to  single  finite  integrals  which  are  more  suitable 
for  numerical  computation.  However,  the  temperature  and  the  thermal 
stress  fields  are  quite  complex  and  remain  in  the  form  of  multiple 
integrals,  for  which  numerical  evaluations  are  obtained  using  a  Simpson 
adaptive  scheme  in  conjunction  with  a  Gaussian-type  quadrature  (Appen¬ 
dix  II).  This  numerical  scheme  is  proven  to  be  effective  since  it 
checks  with  the  numerical  result  for  the  temperature  field  obtained 
from  the  Green's  function  method  from  which  the  resultant  integrals  are 
much  easier  to  evaluate  numerically. 

3.2  Mechanical  Stress  Field 

Basic  Governing  Equations.  The  governing  equations  come  from  Cauchy's 
law  and  Hooke's  law.  In  terms  of  the  moving  convective  coordinates 
{xi }  and  with  the  absence  of  body  forces,  the  acceleration  in 
Cauchy's  law  will  have  only  convective  terms.  They  are: 

3j  ■  pV2  Juat  .  (3.1)* 

4Summation  convention  is  used  for  repeated  indices  of  roman  minuscules. 


I.j.k  *  1.2,3 


(3.2) 


*1j  *  X3kVlj  +  u*3iuj  +  VV  • 


where  Is  the  displacement  field,  Is  the  mechanical  stress  field, 
denotes  partial  derivative  with  respect  to  x^,  5.^  is  the  Kronecker 
delta,  V  Is  the  asperity  traverse  speed  and  p,  x,  and  y  are  material 
constants.  The  method  used  to  solve  the  equations  above  is  similar  to 
the  one  used  by  Eason  [21]. 

Equations  (3.1)  and  (3.2)  are  conveniently  solved  by  the  method  of 
double  Fourier  transform  with  respect  to  (x^  x2) 
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With  the  use  of  the  relation 
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where  the  transformed  quantities  are  denoted  by  a  superposed  (”).  R3 
Is  the  traction  on  the  surface  x3  *  0,  D  denotes  derivative  with 
respect  to  x3  and  the  index  i  »  1,2,3. 


Displacement  and  Stress  Fields.  The  solution  of  the  set  of  ordinary 
differential  equations  (3.4)  and  (3.5)  is  obtained  by  combining 
Equations  (3.4)  and  (3.5)  with  the  result 
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in  which 

M1  *  UT^i)  1/2  :  M2  3  v<^1/2  J  02  -  (X  +  2u)/y  . 

Equation  (3.8)  is  a  system  of  three  simultaneous  homogeneous  linear 


equations  for  three  unknowns  {u^}.  For  nontrivial  solutions  the 
determinant  of  coefficients  should  be  zero  identically,  that  is 
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where  the  distorted  distances  are 
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Hence,  Equation  (3.8)  results  in 
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whose  solutions,  with  the  consideration  of  the  regularity  condition 
(3.7),  are 
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where  the  coefficients  {A^ ,  ,  C^,}  are  independent  of  x3. 


Substituting  Equations  (3.10)  Into  (3.8)  gives  certain  relations 
between  these  quantities  as  follows: 
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A3,  Blt  and  B2  are  to  be  determined  by  the  surface  traction  condi¬ 
tions.  Thus,  Equation  (3.10)  can  be  rewritten  as 
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Expressions  for  the  transformed  stress  components  may  be  obtained 
by  the  substitution  of  Equations  (3.11)  Into  (3.5).  In  particular 
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''.v/.*,  ’’o**  *  **•/  »/  •/  s'  s'  \*  s'  %’A 


In  view  of  the  boundary  condition  (3.6),  A3  and  Bj,  B2  can  be 
readily  solved  and  expressed  In  terms  of  T^*. 


2yFA3  *  n^[(x^  +  7^  -  j  MgX^)  T?3  +  in^x^  +  x^)] 


2uFB1  =  -1x1n1n2^3  +  i-  {[7^(7^  +  x*  -  j  M^x^  -  2n1n2)  +  n2(x*  +  7^ 


-  7  Mg*?)]  \  "  +  7^  -  j  M2x^  -  2n1n2)  IT2}  , 


2uFB2  *  -1x2n1n2R3  +  i-  { “x ix 2 +  x2  “  7  M2X1  ”  2nin2^  ^1 


+  [7^(7*  +  “  7  m|xJ  "  2n1n2)  +  n2(x*  +  7^  -  j  M^)]  JT2}  . 

(3.13) 

Following  the  substitution  of  Equations  (3.13)  Into  (3.11) 

_  ^  o  ^  —  "H  ^ 

2uF  ux  *  (iXjn^  -  x^TTj  -  x^n^)  e  -  Ix^n^ 

"2  m 

x2  xi  Xp  )  ~n?x-> 

-  (n3  •  2nln2>  *  n2n3 J  *1  *  —  (n3  -  2nln2)  R2  j  e  • 


_  —  ^  —  o  i  x«  |  ^ 

2yF  u2  =  1 1  x2  n3T?2  -  7^  -  7|n272)  e  -  (Ix^n^. 


x.x2 

+  -i-i  (n.  -  2n 
n2  j 


ft  1 

ln2)  TTj  “  I  n -  (n3  -  2nxn2)  +  n2n3 


IT,  e 


)  -n^ 


•  ••'.•■‘I  •' N\'\ 


2yF  Oj  a  (ni +  Ix^n^  +  ix^n^)  e 


—2  —2  —  _  ox3 

-  [n^T^  +  x|)  £3  +  iXjH^  +  ix2n372]  e 


(3.14 


where 


2  o  1  0  2  2  2  2 

n3  a  x^  +  x|  -  j  M2  7^  ,  F  =  nj  -  nin2^*l  +  x|)  . 

Stress  components  can  then  be  determined  by  the  substitution  of  (3.14) 
into  (3.5)  in  accordance  with  the  new  variables  as  follows: 


*  n3  —  2n^n2  j  N2  a  h^/n2  |  B  s  ^3*^3  ^2(x^R^  ^  x2R2) 


Thus  we  have 


2UF  0^  -  B  C(x^  +  x2  -  n3)  \  +  2  x^]  e  1  3  -  2  fx^n^ 


+  IXjCx^Ng  +  n2n3l  -  i7^x2N272}  e 


-n2x3 


J4|  ^  Q  o  p  o  1^-3  p  . 

2yF  o22  a  B  C(xJ  +  x2  "  nj)  *  +  2ux2]  e  ■  2u{x2njn2R3 


-  fXjT^N^  +  i72C^N2  +  n2n3]  72}  e 


-n2x3 


2ur  ^33  *  B[(7*  +  x2  -  n*)  \  -  2unJ]  e 


2,  ”nlx3 


+  2„n0[(7?  +  xf)  n.7,  +  in-Jx.lT,  +  xJT.,)]  e 


-n2x3 


2uF  o^2  3  Zyx^B  e  1  J  -  +  jr  x^Cr^n-j 


- 7T  "nlX3 


i  — 


-  (x2  -  x2)  N2]  rx  +  j  x1Cn2n3  +  (x2  -  x2)  N£3  1T2}  e  23 


”nlx3 


2yF  0^3  3  -2i»ix2njTr  e  1  J  +  yfix^O^  +  x^  +  n2)  T?3 


-  2XjX2n^n2lT^  -  [x^  +  n3(n2  +  72)]  T^}  e 


-n2x3 


2yF  3  -2u1x1n1‘B‘  e  1  3  +  yflx^^x^  +  x^  +  n2)  lf3 

-  [x^  +  n3(n2  +  x2)]  ^  -  ^xx2nln2R2^  e  2  3  •  (3.15) 

The  general  solutions  of  displacement  and  stress  fields  are  left  in  the 
Fourier  transformed  expressions  (3.14)  and  (3.15).  The  actual  fields 
are  to  be  obtained  through  an  Inverse  transform  after  the  substitution 
of  the  specific  boundary  values  {TTf } . 


3.3  Temperature  Field 

The  heat  equation  with  constant  thermal  properties,  assuming 
quasi-steady  state  and  no  heat  generation  in  the  medium,  as  expressed 
in  terms  of  the  convective  coordinates  {x-f } ,  is 


3iiT  3  7  31T 


(3.16) 


where  <  is  the  thermal  diffusivity. 
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The  boundary  conditions  at  x3  *  0,  are 


J-qtXj,  x2)  ,  in  the  contact  region 

(3.17) 

0  ,  elsewhere 

T,  3^T  0  as  (x.-Xj)1^2  •  .  (3.18) 

Note  that  T  is  the  temperature  above  the  stress-free  ambient,  k  Is 
the  thermal  conductivity,  q  is  the  heat  flux  (-RjV  =  ufpV),  where  jjf  is 
Coulomb  coefficient  and  p  is  normal  pressure  distribution  on  the  sur¬ 
face  x3  *  0.  It  Is  assumed  that  the  mating  seal  which  contains  the 
asperity  is  an  insulator. 

The  double  Fourier  transform  of  Equations  (3.16)  through  (3.18) 
become,  respectively, 

CD2  -  (s2  +  in)]  T  *  0  ,  (3.19) 


!-<q(7j,  72)  ,  in  the  contact  region 

(3.20) 

0  ,  el sewhere 

T  ♦  0  as  x3  +  -  ,  (3.21) 

where  s'  =  +  x2  and  n  *  -(V/«)  7^.  Mow  let  T  *  ^  +  iT2  and  Q/k  * 

?!  +  i72,  then  Equations  (3.19)  and  (3.20)  result  in 

(02  -  s2)  T1  *  -n?2  , 
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and 


(3 


(D2  -  s2)  T2  3  nT: 


.2 


DT1  “  P1  * 


DT2  3  ?2  a<t  x3  s  0 


(3.2 


Also,  (3.21)  becomes 


Tj,  T2.0  as  x3  ♦  •  . 


(3.2 


Equation  (3.22)  can  be  expressed  as  a  fourth-order  homogeneous 


ordinary  differential  equation  which  is 


[D4  -  2s2D2  +  (n2  +  s4)]  Tp  3  0  ,  r  =  1,2  . 


(3.2 


In  view  of  the  regularity  condition  (3.24),  the  solution  of  Equation 
(3.25)  has  the  form 

-wx. 


Tr  3  e  [Ar  cos(ex3)  +  Br  sin(ex3)] 


(3.2 


where 


^  (s4  +  n 


2)1/2  _  s2 


1/2 


f<s4  +  n2)l/2.S2l1/2 

■  - 2 - 


Substituting  (3.26)  into  (3.22)  and  (3.23)  gives 


P_  =  -uiA_  +  eB_ 


r  *  1,2 


(3. 


0 


1 


CT*2  * 


v  o .  e.'a } 


Hence,  we  have  four  equations  for  four  unknowns  A1#  A2,  Bt  and  B2, 
which  can  be  readily  solved. 


h" 


2a»2'P‘1  +  n?2 

2""  5"  * 

2<i)(ti>  +  0  ) 


n?1  -  2«2P2 

*2  2«(u2  +  0*) 


28  r.  - 

5l - 2 - T“  * 

2e(ai  +  0  ) 


riF^  +  202?2 


2  20U*  +  0*) 


(3 


Following  the  substitution  of  the  expression  (3.30)  Into  (3.26)  and 
according  to  T  *  Tj  +  iT2,  we  obtain  T  as 

_  -wx,  _ 

T  *  e  0  f[C1  cos(0x3)  +  C2  sin(0x3)]  Pj^ 


■■A 


+  [C3  cos(0x3)  +  C4  sin(0x3)]  P2}  , 


(3 


Again,  the  temperature  field  Is  left  In  the  Fourier  transformed  expres¬ 
sion,  which  will  be  used  in  solving  the  thermal  stress  field.  To  be 
complete  and  also  for  comparison,  the  Green's  function  solution 
approach  for  the  temperature  field  is  given  In  the  following. 

The  Green's  function  solution,  following  Carslaw  and  Jaeger  [47], 
can  be  expressed  as 


g(Xj  —  Xj ,  m  x2»  x3 »  t  —  t  )  — 


4pc[iric(t  -  t' )]' 


(  C(x.  -  x|)  -  V ( t  -  t')]2  +  (x2  -  x’)2  +  x2 

•  «P  j  -  — — ■■ - - 1  (3‘32' 

which  is  the  temperature  at  time  t  at  the  location  (xlt  x2,  x3)  due  to 

»  I 

the  unit  heat  flux  input  emitted  at  time  t'  at  the  location  (xlt  x2, 

0).  It  can  be  easily  shown  that  Equation  (3.32)  satisfies  the 
governing  heat  equation  (3.16). 

The  temperature  field  can  be  obtained  by  the  convolution  integral 


T(x^ ,*2 »xj » 3  J*  Kg(x^  “  x^*  x2  ”  x2»  x3»  ^  t  )  dx^  dx2 


where  B^  indicates  the  heat-flux  input  region  and 


/t 

g(x^  -  x^,  x2  -  x2,  x3>  t  -  t')  dt'  . 


(3.33) 


(3.34) 


After  substituting  (3.32)  into  (3.34)  and  changing  the  variable  by 
letting  t  *  R/2/<(t  -  t’ ),  K0  becomes 


VtXj-XjJ^K 


kR  -'R/2/^t 


e-T2-(Y2R2/16ic2T2)  dT 


(3.35) 


where  p  a  density 


c  »  specific  heat 

k  *  pC< 

R  *  [(Xj  -  x')2  +  (x2 


xj)2  ♦ 


For  the  steady-state  solution,  i.e.,  t  Equation  (3.35) 


results  in 


G(x,  -  xi,  x2  -  x2,  x3)  =  lim  Kq 

A  t-H* 


l  -V[R  -  (xx  -  x^)]/2< 

OTTe 


(3.36) 


v,  A  **.  •  .  •*. 

V  V  V  -  *.•'  *.•  «  v.  *  -  •  .  •  *  *  .  »  .  v*.»  *  .  ■  «  .  •  . 


Accordingly,  the  steady-state  temperature  field  can  be  obtained  as 


T<WX3>  -s  L  q<Xl,x2) 


-VCR-(Xl-x^)3/2< 


dx' 


dx^ 


(3.37) 


which  will  be  employed  to  give  numerical  results  for  different  cases  of 
application. 


3.4  Thermal  Stress  Field 

The  governing  equations  for  the  thermal  stress  solutions  are  based 
on  the  theory  of  quasi-static  uncoupled  thermoelasticity.  The 
governing  equations  are  then 


U^Uj  +  (x  +  u)  3^(3,^)  *  (3*  +  2y)  a3jT 

(3.38) 

and 

“ij  *  xWtj  *  +  W  ■  <3x  +  2“>  “Tsij 

(3.39) 

subject  to  the  boundary  conditions 

«31  "  0  *  1  ’  1*2*3  at  x3  *  0  ’ 

(3.40) 

T  1/2 

aij  »  uj  ■”  0  as  (xjxi^  •  • 

(3.41) 

The  double  Fourier  transform  of  Equations  (3.38)  through  (3.41)  are 


-  IagXjOUj  -  ia^DUg  +  U^Q  -  s  )  u3  *  a3DT  , 


——  »  -iUjlTjU^  +  a^I Tg)  +  a^Du^  -  a^T  , 


— jp  *  -i(a4x^  +  a ^ITg)  +  a^Du^  -  a3T  , 


—  »  -Ta4(x1u1  +  72u2)  +  aiOu3  -  a3T  . 


— - 1(*2“l  *  xl“2>  • 


-T 

°23  — 

—  a  Qu- 
u  2 


^x2U3  ’ 


-T 

°31  — 

—  *  Du. 
u  1 


ixlu3  » 


"Ii  ■  0  • 


1  »  1,2,3  at  X3  *  0  , 


— T  —  —  —  1/2 

o<  i  ,  u.  ♦  0  as  (x<x. ) 


where  a^  *  x/u,  al  ■  a4  +  2,  a2  a  a4  +  1  and  a3  ■  (3au  +  2)  a. 

Having  a  known  temperature  expression  as  In  Equation  (3.31),  we 


may  choose  the  solution  for  the  displacement  components  I?!  and  U2  in 


the  form 


*  f(x3)  e  , 


U2  a  g(x3)  e 


(3.53) 


where  functions  f  and  g  are  to  be  determined  later  by  combining  Equa¬ 
tion  (3.43)  and  Equation  (3.42)  which  results  in 

(D2  -  s2)(72u1  -  x’jITg)  *  0  . 


Considering  Equation  (3.53),  it  follows  that 


f  a:r  9  . 


Equation  (3.42)  can  be  rewritten  as 


(3.54) 


Du 3  *  — — -  C-ia3x1T  +  (a.^  +  x2)  -  D  u^  +  a^x^] 


(3.55) 


In  view  of  D2^  obtained  by  taking  the  differentiation  of  Equation 
(3.55)  with  respect  to  x3,  it  follows  from  Equation  (3.44)  that 

u3  a  — Z — J  C -a ilT^  +  UjS  +  a2xp  Du^  +  a2X]X2Du2  *  1a3x .(3.56) 
a2xls 


Differentiating  u3  with  respect  to  x3  and  then  equating  the  result  with 
Equation  (3.55),  the  00E  governing  the  displacement  field  becomes 

-a jD4^  +  (a2x2  +  s2  +  a^2  -  a2x2)  D2^  -  s2(a1x2  +  x2) 


+  agXjXgD2^  -  a^x^x^  =  ia-jX^D2  -  s2)  T  .  (3.57) 


/*\yV*v*I/-V'y*y  *  *  y-y*y‘\-y\V 
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From  Equations  (3.31),  (3.53),  and  (3.54),  Equation  (3.57)  leads  to 
-D4g  +  4<uD^g  +  2{s2  -  3u)2)  D2g  +  4tu02Dg  -  s4g 

2a3  - 

*  aj“  x2weC(B2  -  18^  cos ( 0  x3)  +  (-A2  +  iA1)  sin(a  x3)]  (3.58) 

where  Alt  A2,  Blf  and  82  are  the  same  as  in  Equation  (3.30). 

The  complementary  solution  of  Equation  (3.58)  (gH)  is  obtained 
by  solving  the  corresponding  characteristic  equation  which  is 

-m4  +  4ojm^  +  2(s2  -  3w2)  m2  +  4ti»82m  +  84  *  0 
or 

(m  -  «  -  s)2  (m  -  «  +  s)2  *  0  .  (3.59) 

Taking  account  of  the  regularity  condition  the  complementary  solution 
is 

u  (urs)  x« 

gH  a  (hx  +  h2x3)  e  (3.60) 

where  hx  and  h2  are  contants  of  Integration. 

The  particular  solution  of  Equation  (3.58)  (gp)  is  determined  by 
the  method  of  undetermined  coefficients  as  follows. 

g  ■  gy— 9  CG1  cos( 8  x3)  +  H1  sin(e  x3)]  (3.61) 

where  Gi  *  B2  -  IB!  and  Hj  *  -A2  +  1At. 

With  gH  and  gp  known,  the  particular  and  homogeneous  solutions 


of  uj  can  be  readily  obtained 


uj  *  ^7  (hx  +  h2x3 )  e 
x2 

-c  "sx3 

*  (hj  +  h2x3)  e  , 

.  -sx- 

»  — L_  [a^O^  +  h2x3)  +  (a2  +  2)  h23  e 

a2X2 


_n  3,X„  “0>X-» 

uj  *  [Gr  cos(e  x3)  +  Hr  sin(e  x3)]  e  ,  r  *  1,2  , 


“S  ■  -tqzs  CG3  cos(9  x3>  +  H3  s1"<9  x3)]  9 


-0»x- 


where 


a  G2  ,  =  ^2 


G3  ■  -iUGj  -  9H1)  ,  H3  *  -KmHj  +  9Gx  ) 


As  a  result,  the  general  solutions  are 


■■  *“0  * C 

ui  3  “i  +  u1 


(3.62 


(3.63 


The  two  unknown  constants  hx  and  h2  will  be  determined  by  substituting 
Equation  (3.63)  into  Equations  (3.47)  and  (3.49)  and  employing  the 
boundary  condition  (3.51).  Thus  we  have 


—  Ca1s9  H1  +  (s  -  +  a1et)  G^]  , 


28^2$  0)6 


&  X 

h2  *  2a^sle  ^"s0  H1  “  s“  +  al8  *  Gl^  *  *3, 

Following  the  substitution  of  the  displacement  solution,  (3.62) 
through  (3.64),  into  Equations  (3.45)  through  (3.50),  the  thermal 
stress  field  can  be  readily  obtained. 

flii  f  ”ti)X^  “sx3 1 

u  *  1  ('blHl  +  b2Gl}  e  +  (‘b3Hl  +  b4Gl)  e  • 

Oyy  f  -SX-1 

-jp-  =  1  (-b5H1  +  bgGj)  e  J  +  (-byH1  +  bgGj)  e  J  , 

r  “six.  ”sx. 

-f  -  1^(-b9H1  +  bloGl)  e  +  (buHl  -  bl2G1)  e  3j  , 

«T,  r  “«■#,  "*x,-i 

~  *  iXlX2 ^”b13Hl  +  b14GlJ  e  +  ("b15Hl  +  b16Gl}  e 


-T  r 

-f -  -T4(t>i7Hi '  biaGi 


”ti)X«  ”SX^ 

1  e  -  <b19Hl  -  “so6!1  e 


—  •  -x,  (b,,H,  -  b,aG,)  e  J  -  (b,ftH,  -  b„nG, )  e  ,  (3. 


3  —2 

—  [2u0  cos(ax3)  -  x1  sin(ax3)]  , 

a3  o 

j^g-  Cxj  cos(ex3)  +  2u)9  sin(ex3)]  , 
\a2s  aia2  ax  /» 


~  [2(00  cos(gx3)  -  x2  s1n( ex3 ) ] 


*3  g 

—  [x£  cos(ex,)  +  2(d9  sin(9x,)]  , 


i  (!L  + 

\  a2s  ala2  al  /  ’ 


*3  rai«  -  S  _2  r(S(D  -  s2)  .  s  rrZ 

^9  a^s  ’  x2  + - ^ - +  (s  u)  x2 


si n ( 9xn)  , 


COS(9X3) 


Note  that  Pj  and  P2  are  the  real  and  Imaginary  part  of  Q/k. 
Therefore 


?!  3  2^-  J  Q(xl#  x2)  cos  (XjXj  +  72x2)  dXj  dx2  , 


F2  Q(X1>  x2>  Sln  (X1X1  *  V2>  dxl 


dx. 


(3.6 


H 


For  the  case  P^  3  0,  Equations  (3.65)  become 
-T 


^■Ti[{bici 


-sx, 

+  b2C2)  e  J  +  (b3C1  +  b4C2)  e  J 


]■ 


-T 


'-f  *  Ti[(bsci 


-04X-.  -sx, 

+  b6C2)  e  4  +  (b7C1  +  baC2)  e  J 


]■ 


-T 


'33 


ri[(b9Cl 


-wx,  -sx- 

*  b10C2)  e  -  (buCx  *  buC2)  e  J 


]• 


— T 


“ir  *  ^ixix2  _(bi3ci +  buc 


-tux-  -sx, 

2)  e  +  +  t»lgC2)  e 


]• 


-T 


*T*  •<x27l[(1)I7Cl  +  b18c 


-ux,  -sx- 

2)  •  -  +  b20C2>  e 


]• 


-T 


V  *  -,xl1S'1[(b17Cl 


-tuX7  -SX- 

+  b18C2)  e  J  -  (b19C1  +  b20C2)  e  U  .  (3.6' 


1  .  (3.( 


3.5  Numerical  Results 


Previous  sections  have  shown  the  general  expressions  for  the 
transformed  mechanical  and  thermal  stress  fields.  After  applying  the 


Inverse  Fourier  transform  and  following  through  the  simplification  of 
the  multiple  integrals,  results  for  several  cases  can  be  obtained 
numerically. 


3.5.1  Rectangular  Contact  Area  with  Uniform  Pressure  (Case  1) 


Mechanical  Stress  Field 


The  tractions  of  the  mechanical  loading  are 

*32  3  “R2  a  0  »  *31  ”  "l"!’  *2* 


*  -R^x^,  X£i  0)  i  a  8  0  »  i  x^»  0) 


where 


pQ  =*  constant  , 


Xj^  £  a,  x2  <_  b 


0  , 


el  sewhere 


R1  *  "wfR3  * 


The  corresponding  transformed  boundary  conditions  can  be  readily 
obtained 

KTjXj+IgXg) 


TT3(xlt  x2,  0 f 

— o»  ^  •• 


R3  e 


dXj^  dx2 


,  /»b  i  (7,  x.  +X,x, ) 

hL  L? o * 


/•b  ix2x2 

J  .  e  d*2 

/•a  ix.x. 

J  e  dxl 

_ J  -b 

.  "a 

2?. 


sin(x,a)  si n Cx-b)  . 


By  following  the  substitution  of  (3.68)  into  (3.15)  and  applying  the 
inverse  transform,  the  mechanical  stress  components  can  be  expressed  as 
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second  the  frictional  loading  T^.  The  expressions  are 
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Equations  (3.69)  can  be  reduced  to  single  integrals  by  using  the 
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The  term  JQ  —  e  ds  cancels  out  during  simplification.  Thus, 
stress  fields  are  expressed  as 
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The  integrals  appearing  in  Equations  (3.70)  are  suitable  for 
numerical  evaluation  at  any  point  Inside  the  solid.  However,  most 
integrals  are  not  defined  at  4  *  0,  */2  or  ir.  Limiting  processes  are 
Imperative  so  that  singularities  can  be  removed. 


Temperature  Field 

Replacing  Q(xltx2)  in  (3.66)  by  expressions  (3.20)  and  (3.17),  the 
heat  flux  input  In  the  transformed  space  is 
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Expressions  above  are  substituted  Into  Equation  (3.31)  for  the 
transformed  temperature  field.  The  inverse  Fourier  transform,  after 
certain  slmplications,  results  in: 
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Alternately,  the  temperature  distribution  can  also  be  obtained  by 


Green's  function  method  which  will  serve  as  a  double  check  for 
numerical  results.  The  temperature  field  results  by  substituting 
Equation  (3.35)  into  (3.33)  and  letting  t  ♦ 
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This  integral  expression  is  more  conducive  to  numerical  computation. 


Thermal  Stress  Field 

Following  the  substitution  of  the  expression  into  Equations 
(3.67)  and  applying  the  Inverse  Fourier  transform,  the  thermal  stress 


field  is  thus 
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Due  to  the  complexity  of  the  expressions,  no  existing  automatic  inte¬ 
gration  scheme  for  the  integrals  above  is  efficient.  Hence  the 
Gaussian-Laguerre  quadrature  formulas  and  Simpson's  integration  scheme 
are  used  for  the  numerical  approximation  (Appendix  II). 


3.5.2  Disk  Contact  Area  with  Uniform  Pressure  (Case  2) 
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The  transformed  boundary  conditions  can  be  obtained  by  using  Equation 
(3.3). 
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Following  the  same  procedure  as  in  Case  1,  expressions  are 


obtained  for  the  stress  components  of  the  forms 
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The  Integrals  for  St  and  Wt  can  be  solved  analytically  so  that  the 
stress  expressions  have  the  form  of  single  Integrals  which  can  be  eval¬ 
uated  numerically. 


Temperature  Field 

?!  and  P‘2  are  given  by 
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Substitution  of  Tl  and  J2  ^nt0  Equation  (3.31)  and  again  applica¬ 
tion  of  the  inverse  transform  results  in 
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for  which  the  numerical  evaluation  can  be  done  by  double  Simpson's 
schemes.  Since  Equation  (3.76)  Is  Integrated  within  finite  Intervals 
It  Is  advantageous  to  compute  the  numerical  results  through  Equation 
(3.76)  Instead  of  (3.75)  where  the  Integral  Is  Improper. 


Thermal  Stress  Field 

In  the  same  fashion  as  In  Case  1,  It  can  be  shown  that 
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3.5.3  Disk  Contact  Area  with  Paraboloidal  Distribution  of  Pressure 
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The  expressions  of  the  stress  components  are  similar  to  those 
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(3.74)  in  Case  2  except  St  and  Wt  being  replaced  by  St  and  Wt  which  are 
of  the  forms 
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Again  S  and  W  can  be  solved  analytically  and  accordingly  numerl- 
cal  results  of  the  stress  components  can  be  readily  obtained  by 
automatic  Integration  schemes  for  single  Integrals. 


Temperature  Field 
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The  temperature  field  can  be  obtained  by  replacing  f2  In  Equations 
(3.75)  by  f3  which  has  the  following  form 
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Thermal  Stress  Field 

The  expressions  have  the  same  form  as  Equation  (3.77)  In  Case  2 


except  f2  being  replaced  by  f3. 


As  an  Illustration  of  the  numerical  evaluation  of  the  Integrals, 
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The  temperature  rise  from  the  cold  state  is  4  a  Tk/qQa,  where  qQ  Is  the 
heat  flux  due  to  the  average  frictional  load,  that  Is,  qQ  *  u^0v* 
Figures  3.2  to  3.4  show  the  mechanical  principal  stresses  (oj!)  on  the 
trailing  side  with  different  depth,  c,  for  Cases  1,  2,  and  3,  respec¬ 
tively.  It  Is  noticed  that  stresses  rise  sharply  in  the  neighborhood 
of  c  *  1,  as  c  approaches  zero  for  Cases  1  and  2.  Such  stress  singu¬ 
larities  at  5  *  ±1  are  expected  because  of  the  discontinuities  In  the 
loading  functions.  Temperature  fields  for  different  cases  are  shown  In 
Figures  3.5  to  3.7,  In  which  the  temperature  decays  rapidly  with 
respect  to  the  depth.  Consequently,  thermal  stresses  are  compressive 
only  in  a  very  thin  surface  layer  due  to  the  thermal  effect.  The  mag¬ 
nitude  of  compressive  stresses  also  decreases  rapidly  with  the  Increase 
of  the  depth  until  it  becomes  tensile.  Figures  3.8  to  3.10  show  the 
corresponding  thermal  principal  stress  (<jj)  at  depths  c  *  10" 1  and 
1.0,  where  tensile  stresses  dominate,  for  each  case.  The  combined 
effects  from  mechanical  and  thermal  fields  for  three  cases  are  given  in 
Figures  3.11  to  3.13.  Accordingly,  the  critical  asperity  pressure  to 
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Maximum  mechanical  principal  stresses  In  the  surface  layer  (Case  2 


Maximum  mechanical  principal  stresses  In  the  surface 


Maximum  thermal  principal  stresses 


Figure  3.12.  Maximum  combined  thermomechanical  principal  stress  field  (Case  2). 


Figure  3.13.  Maximum  combined  thermomechanical  principal  stress  field  (Case  3). 


Initiate  cracks  can  be  predicted.  For  an  asperity  with  a  size  of 
0.5  mm  (0.02  in.),  a  velocity  of  15  m/s  (600  ips)  and  a  coefficient  of 
friction  »  0.5,  the  maximum  principal  stress  for  Cases  1,  2,  and  3 
reaches  the  ultimate  strengh  at  critical  pressures  of  353  MPa  (51,253 
psi),  318  MPa  (46,083  psl),  196  MPa  (28,400  psi),  respectively.  In  al 
cases,  heat  checking  by  the  mechanisms  postulated  may  occur  on  the 
trailing  side  at  a  depth  of  0.025  irm  (0.001  in.)  from  the  surface. 


CHAPTER  4 


A  THREE-DIMENSIONAL  MODEL  OF  PERIODIC  MOVING  ASPERITIES 

4.1  Mathematical  Model 

The  established  three-dimensional  theory  of  a  single  moving  asper¬ 
ity  in  Chapter  3.0  is  expanded  to  incorporate  the  more  realistic  check¬ 
ing  failure  in  a  friction  loaded  medium  shown  In  Figure  4.1.  As 
demonstrated  by  Burton  [7],  the  hot  spot  at  the  Interface  between  a 
stationary  and  a  rotating  contacting  surface  precesses  at  a  regular 
rate.  With  reference  to  a  region  on  the  surface,  the  hot  spot  is 
indeed  moving  over  the  surface  at  a  periodic  regularity.  Under  a  sin¬ 
gle  severe  asperity,  the  seal  may  initiate  a  crack.  With  repeated  mov¬ 
ing  friction  loads,  the  crack  may  be  initiated  at  a  much  lower  friction 
load  and  may  be  associated  with  a  different  fracture  configuration. 

The  geometries  of  the  contact  area  and  loading  profiles  under  consider¬ 
ation  will  be  classified  Into  Cases  1,  2,  and  3,  which  are  the  same  as 
those  cases  In  Chapter  3.  Moving  coordinates  fixed  to  the  center  of 
the  first  moving  asperity  and  a  time  coordinate,  (x x  -  x2  -  x3,  t),  are 
used  throughout  the  analysis.  It  Is  assumed  that  the  repeated  asperi¬ 
ties  occur  with  a  period  of  t  and  move  with  the  same  velocity  (Y). 
Therefore,  the  center  distance  between  asperities  ( x ’  *  Vr/a)  will 
later  play  an  important  role  in  illustrating  the  cumulative  effects. 

The  period  t  is  an  Important  parameter  because  of  the  following 
arguments: 

I.  If  the  period  is  so  large  that  the  thermomechanical  effect 
resulting  from  the  first  asperity  has  already  died  out  at  the  time  the 


second  asperity  appears,  then  the  physical  and  mathematical  models  will 
be  those  of  the  single  moving  asperity. 

II.  If  the  period  Is  shortened  so  that  the  mechanical  effect  of 
the  first  asperity  has  dissipated,  then  the  thermal  effect,  which  dif¬ 
fuses  at  a  slower  rate,  remains  to  affect  the  thermomechanical  response 
to  the  following  asperity. 

III.  If  the  period  is  very  short  (high  frequency),  then  both  the 
mechanical  and  the  thermal  responses  of  the  first  asperity  could  Inter¬ 
act  with  the  following  one. 

The  results  In  Chapter  3.0  are  therefore  applicable  to  the  long 
period  excitations.  Category  III  applies  to  loads  traveling  at  wave 
speed.  For  the  problem  under  consideration,  such  speeds  do  appear  In 
realistic  situations.  Consequently,  It  Is  Important  to  consider  Cate¬ 
gory  II  for  the  present  problem.  It  Is  significant  to  investigate 
first  whether  a  bounded  steady  state  solution  exists.  If  it  does,  the 
thermal  stress  field  after  long  term  operation  would  also  reach  a 
steady  state  solution  which  will  be  obtained  analytically  by  Fourier 
transform.  General  solutions  of  different  cases  for  the  mechanical, 
temperature,  and  thermal  stress  fields  will  be  given  In  Sections  2  to  4 
and  numerical  results  will  be  shown  In  Section  5.  The  variables  used 
throughout  this  chapter  will  be  the  same  as  those  In  previous  chapters 
unless  specified  otherwise. 


4.2  Mechanical  Stress  Field 


approach  for  the  mechanical  stresses  are  exactly  the  same  as  those  of  a 
single  asperity  in  Section  3.2.  However,  it  should  be  mentioned  that 
the  boundary  conditions  for  the  current  periodic  loading  problem  are 
modified  to  be 

”31  ■  '  Ri  “  x3  “  0  • 
and 

ui  »  j  ♦  0  as  (x^)1^2  ♦  •  (4.1) 

where 

P(xltx2)  ,  8^:  contact  regions,  j  *  1,...,N 

0  ,  elsewhere 


R1  3  ,  R2  - 


where  P  Is  the  normal  pressure  distribution  prescribed  over  the  contact 
regions,  which  is  quite  general  at  th*s  stage.  N  Is  the  total  number 
of  existing  moving  asperities. 

Following  the  same  approach  and  using  the  same  notation,  the  solu¬ 
tions  subject  to  the  boundary  condition  (4.1)  are  obtained  as  follows. 


2uF^i  =  U[(s2  -  n2)  x  +  2U  x2]  e’"1*3 


-  ZuCx^njn^  +  Ix^x^  +  n2n3^  ^  e 


”n2x3 


-n,  x. 


2UF«22  =  T5C<S2  -  X  *  2UX^]  .  13 


■  2u[x?n,n2^3  *  ff.AjM  e 


"n2x3 


.1 


2yF  0^2  a  Ot(s^  -  n^)  \  -  2y  n^]  e 


2,  "nlx3 


o  -n9x^ 

+  2u  n 2 ( s +  in^XjTT^]  e  0 


?  c  -M 
2yF  a12 


— _  “nlx3  _  i  _ 

2y  Xj^2  '5'  e  -  211(3^  Xjn^TTj  +  j  XgCngn- 


-  (x*  -  x2)  N2]  Rl}  e 


-n2x3 


“nlx3 


2uF  0^2  *  “^yi  X2nl  A  J  +  u{1x2n^(s^  +  n^) 


-  ^xix2nin2^1^  e 


-n2x3 


-  .  ~n  ~nlX3 


2 u F  o^j  *  -2pi  x^  D  e  1  0  +  y{ix^n^($^  +  n 2)  1T3 


-  [x^N^  +  n3(n2  +  x*)]  TTj}  e 


'n2x3 


(4.2) 


where  D  ■  n3l?3  +  1n2xllT1,  the  remaining  coefficients  are  defined  as  in 
Chapter  3.0. 


4.3  Temperature  Field 

It  is  assumed  that  the  temperature  field  has  a  slow  dissipation 
rate  and  remains  at  an  appreciable  magnitude  during  the  passage  of  the 
second  moving  asperity.  In  other  words,  at  every  traverse  of  an  asper¬ 
ity,  the  corresponding  initial  temperature  is  modified  after  the  pas¬ 
sage  of  the  preceding  asperity.  The  analysis  therefore  proceeds  with  a 
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time  sequence,  at  the  period  of  the  asperity  occurrence  (t),  to  solve 
the  transient  temperature  field  of  each  asperity  In  accordance  with  the 


Initial  condition  resulting  from  the  previous  response.  Let  Ta 
denote  the  transient  solution  due  to  the  <*th  moving  asperity  appearing 
at  the  instant  t  *  (a  -  1)  t.  The  governing  equation  for  unsteady  heat 
condition  in  terms  of  the  coordinate  system  (xj  -  x2  -  x2,  t)  has  the 
form 

3„T  *  i  U.T  +  Y3.T  }  ,  t  >  (a  -  1)  t  .  (4.3) 

11  a  <  tala  — 

The  solution  of  Equation  (4.3)  must  satisfy  the  following  Initial 
and  boundary  conditions 

T  *  U  when  t  *  (a  -  1)  r  (4.4) 

<t  CL 


where  Ua  is  the  Initial  temperature  resulting  from  the  preceding 
asperities. 


ka-T 
3  a 


iPtx^.Xg),  Bj!  contact  regions,  j 

-ufV 

(0  ,  elsewhere 


1,2 


(4.5) 


At  x3  *  0, 

T  ♦  0  as  x,  *  • 
a  3 

where  a  3  1,2,...,N. 

Equations  (4.3)  to  (4.5)  can  be  solved  by  Green's  function 
approach.  First,  consider  a  3  1  and  denote  Tj  as  the  temperature  above 
the  traction-free  ambient.  Then,  setting 


Tj  3  u»2  exPt'j^  (*i  "  "j  t)] 
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\*  %*  ■.* 


V  \  N 


we  have 


3ii“i  =  <  V: 


mj  *  0  at  t  a  0  , 


k33“l  *  Qi  *  exPC_  (*i  “  j  t)] 


for  which  the  Green's  function  solution  can  expressed  as 


I  I  I 


G (x^ tXg »t  ) 
and  hence , 


(  C(xl  -  Xl)2  Mx2  -  x;>2  +  x2]| 
4Pct.«(t  Pi  — : — j 


“l  ■  BfVJ^  f  pt*i.xa)  •  exp[-  (Xj  t')]  G(xj,x2,t)  dx]  dx]  dt 


and 


Tl  ’  “fV  fo  Is.  P<Xl,X2>  '  9(Xl  '  xi’ 


•  III 


*2  "  Xg.Xj.t  •  t  )  dx ^  dx 2  d t 


(4.6] 


where 


/  1  l  (xi  ’  Vt)2  +  x2  +  x2) 

exp| - cr-2— 2!  <«•« 


which  is  essentially  the  Green's  function  solution  for  Equations  (4.3) 
and  (4.5). 

Now,  considering  a  *  2,  we  have 


3^2  *  ^  +  V 3 jT 2 ^  *  t  >  t 


K 


Tg  *  U2  at  t  *  t 


(4.8) 


k33T2  3  Q2  3  "wfV  ) 


P(x,,x,)  ,  B,  and  B. 


,  el sewhere 


at  *  0  . 


Mote  that  the  Initial  temperature  can  be  written  as  U2  ■  T1(x1,x2,x3, 
t),  where  TL  is  given  in  Equation  (4.6). 

The  solution  of  Equation  (4.8)  can  be  readily  determined,  which  is 

T2  3  pC  Vxi»x2*x3*t)  *  g(xl  "  X1  ,X2  ’  X2’X3  '  X3,t} 


2  „t 


dx3  dx2  dxi  +  “fV  2  jf  /B  P<xi’x2 


g(x1  -  xx,x2  -  x2,x3,t  -  t  ) 


»  t  I 

dx^  dx2  dt  . 


(4.9) 


Similarly,  the  recurrence  formula  can  be  obtained  for  subsequent  tem¬ 
perature  fields.  Therefore,  the  transient  solution  for  the  time  after 
(N  -  1)  t  has  the  following  form 


TN-’c/  f  TN-1^X1,X2,X3’  ^  "  !)  t)  • 


-  x>  dx;  dX;  *  „fv  i  £1)t/(. 


•  “  *2**3 


1  111 


•  P(x,,x0)  •  g(x,  -  x .  ,x0  -  x _ ,x , t  -  t  )  dx,  dx-  dt 


(4.10) 


We  observe  that  the  first  term  of  Equation  (4.10)  converges  uniformly 
to  zero  as  t  approaches  Infinity.  Hence,  for  the  steady  state  solu¬ 
tion,  Equation  (4.10)  yields 


r  ^  /*t  /•  i  i 

ufy  lim  2  /  /  P(x.,x5)  •  g(x, 

f  t~  U-i  Am-Dt  jb.  1  2  ] 

v 

'  '  '1 

•  dx^  dx2  dt 

Making  use  of  the  succeeding  result 

/•t  it  ii 

lim  /  g(x.  -  x^x,  -  x9,x,,t  -  t  )  dt 
t-^H-Dt  1  12  23 


Xi»x2  —  x2,x3,t 


- 1 ) 


(4.11) 


V  (x  i~x  ^) /2»c 
ir3/^  kR 


lim 

t-*» 


/ 


exp 


R/2/t""(¥-l)T 


H 


y2  ♦ 


v2r2 


16  k  y 


u) 


dy 


*  27kff  exp{~  V^R  ’  (xl  ”  x^)]/2k}  , 
in  Equation  (4.11),  we  can  rewrite  the  steady  state  solution  as 

ufV  N  -  ,  ,  exp{-  V[R  -  (x.  -  x. )]/2k}  ,  , 

T  "  SwF  ?,  L  p(xi»x2)  *  R  dxl  dx2  *  (4,12) 
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The  expression  above,  while  leading  to  the  solution  of  the  temperature 
f1ei«,  cannot  be  applied  readily  in  solving  the  thermal  stress  field. 

An  alternative  method  with  double  Fourier  transform  is  then  used.  We 
start  by  taking  the  transform  of  the  heat  equation,  using  the  corres¬ 
ponding  initial  and  boundary  conditions  for  Tfj,  and  then  simplifying 
the  result  by  the  relation 


results  in: 


9i  »t)  3 


1 


Equations  (4.13),  (4.15),  and  (4.16)  lead  to 


T  *  fi££ 

‘N 


-ie(s2+in)  t  /•• 


(met)172  •'O 


£  ^ (x^ ,X£ ,X3 ,  ( N  “  1) 


:)  •  e 


(x3-xj)2 


UN(xl,x2)  /*t  exp{-[ic{s2  +  1n)(t  -  t  )  +  x2/4<(t  -  t  )]} 


pc(iric)  '  ■'(N-I)t 


(t- t^r 


dt’  . 
(4-17) 


Introduce  the  change  of  variable  as  follows 


z«(t-t)1/2 


then  transform  the  last  term  of  Equation  (4.17)  Into  the  form 

2TLCX.  ,77)  /-Ct-tM-Dx]1^2  -[«(s2+in)z2+x2/4<z2] 

-  — n  -  /  e  3  dz  .  (4.18) 

pcU<)1/Z  J0 


By  letting  t  the  first  term  of  Equation  (4.17)  vanishes.  In 
connection  with  the  expression  (4.18),  it  follows  that  the  steady  state 
solution  can  be  written  as 


-2Hm(x1,x2)  -U(s2+1n)z2+x2/4KZ2] 


Note  that 


(4.19) 


f  xv_1  exp(-  bxu  -  yx"u)  dx  »  ~  (-J)^  K  (2/e7) 

•'0  u  B 

u 

where  K  is  the  modified  Bessel  function  of  the  second  kind  of  order 
u 

Therefore,  Equation  (4.19)  becomes 


u 


Qn  “(ox-j 

T  *  e  lci  cos(0  x3)  +  c2  S1’n(0  x3^  (4.20) 

which  is  the  steady  state  temperature  field  due  to  N  periodic  moving 
asperities  and  will  be  used  to  derive  the  corresponding  thermal  stress 
field  in  the  following  section. 


4.4  Thermal  Stress  Field 

The  governing  equations  and  boundary  conditions  are  the  same  as 
given  in  Equations  (3.38)  to  (3.41),  their  transformed  expressions  are 
also  the  same  as  Equations  (3.42)  to  (3.52).  In  view  of  the  tempera¬ 
ture  expression  (4.20)  which  Is  similar  to  the  temperature  Equation 
(3.31)  for  a  single  moving  asperity,  the  solution  approach  in  Section 
3.2  can  be  applied.  However,  a  different  method  of  Inverse  Differen¬ 
tial  Operator  is  presented  here  so  that  we  can  check  the  solutions  from 
both  methods. 

Equations  (3.42)  through  (3.44)  form  a  system  of  three  simultan¬ 
eous  linear  differential  equations.  The  dependent  variables  Hi  are 
to  be  solved.  The  system  can  be  rewritten  as 
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Ali^  * 

Ai  , 

i  *  1,2,3 

where  a  is 

the  determinant 

of  the  system. 

D2  -  (a^  +  x2) 

“a2*l*2 

-iagXj^D 

A  *  j 

-a  2X^2 

D2  - 

*  (ajx|  *  *i> 

•iagJg0 

-iag^D 

-iagXgD  a 

jD2  -  s' 

*  < 

a2(D2  -  s2)3 

and  Ai's  are 

-1a3XiT 

-*2*1*2 

-**2*1° 

A1  * 

-ia-^T  D2  - 

<al*2  * 

—2  — 
x^ )  -ia2x2D 

a3DT 

-ia272D 

./  -  s2 

-ia,7.  (D2  -  s2)2 

r , 

cation  of  Cramer's  rule 


Rewrite  Equation  (4.20)  and  Its  derivative  with  respect  to  x3  as 


-UX3  1 


T  *  e  [F  co$(ex3)  +  G‘  si n{ ex3 ) 3 


and 


“04X3  1  1 

DT  »  e  lF^  cos( 0  x3)  +  G^  s1n(e  x3 )3 


(4.2: 


respectively,  where 


•  "(Tu 


.  Q, 


F  *  r  ci  • 


T  u2 


c~  , 


F^  *  8G  -  wF  ,  *  “(uG  +  0^)  • 


In  accordance  with  Equation  (4.23),  we  solve  Equation  (4.22)  as  fol¬ 
lows. 


uj  «  — .  (D2  -s2)’1  |e  3  [F  cos(ex3)  +■  G  sin(e  x3)]| 


-1a,X,  ”UX-j  n  J  .1  •  * 

— r=-=>  e  *  •  (D*  -  2uD  +  8*)  1  [F  cos(0x3)  +  G  s1n(e  x3)] 


‘1 


1a«jX.  ^  1  1  1 

Jie  J  •  D  1  [F  cos(0  x,)  +  G  sin(8  x.)] 


2a 


spr ca'  cosU  x3>  * F'  s,n(9  *3)] 


(4.2 


Note  that  use  has  been  made  of  the  fact  that 


2  .  e2  .  a2 

<i»  *  S  +  8  . 


In  a  similar  fashion  we  obtain 


“2  ■  Tl i^T  CG'  “s(9  x3>  '  F'  s1n(9  x3^  • 


p  o  p  .i  l  l 

“3  *17  •  (D  ’  S  )  e  tFl  C0S(9  X3}  +  G1  Sin(9  x3)] 


^3  1  1 

CG1  cos ( 8X3 )  -  Fx  si n(  9X3)3  . 


(4.25) 


(4.26) 


The  complementary  solution  is  determined  by  au-,-  =  0  which  is 


equivalent  to 


(D2  -  s2)3  (U^.L^.U^)  *  0 


of  which  the  solutions,  with  the  consideration  of  the  regularity  condi¬ 


tion,  are 


uj  »  (aJ  +  B^.x3)  e  SX3  ,  i  =  1,2,3 


(4.27) 


where  A/s  are  unknown  constants  which  are  to  be  determined  by  boundary 

I  I 

conditions.  B^'s  are  expressed  in  terms  of  A/s  as  shown  below. 


■  a?  x.  1  1  1 

B.  *  -  - — — — *  •  --  (7. A.  +  x_A0  -  isA,) 
l  a2  +  2  s  11  22  3 


32  “  *  a2  +  2  *  T  (xlAl  +  x2A2  “  1SV  * 


B3  a  1  a2  +  2  ^X1A1  +  x2A2  "  1  sA3 ^ 


(4.28) 


According  to  Equations  (4.24)  to  (4.27),  the  general  solutions, 

_  _r  _ p 

u^  ■  uj  +  ulj,  have  the  form 

i  i  *SX«  i  i 

^  *  (A^  +  B..x3)  e  +  [Ci  cos(e  x3)  +  0^  sin(e  x3)]  e  (4.29) 
In  which 


.  -ia.,x  1 

cr’2q3rG  • 

•  la-sX  < 

Dr-IifsF  •  r*1-2 

1  4.  1 

1  ”81  1 

C3  *  G1  » 

d3  Iaju0  F1 

In  view  of  Equations  (3.47),  (3.49),  and  (3.50)  and  taking  account 

traction-free  boundary  conditions  leads  to  a  system  of  three  slmul- 

taneous  equations 

dHj  -  17^3  »  o  , 

Dug  “  *  0  , 

-ia^^u^  +  72^2?  +  a]_D73  -  a3T  »  0  .  (4.30) 

Substituting  expression  of  T  (4.23)  and  "u-j  (4.28)  into  Equations 

I 

(4.29)  to  solve  three  unknowns  A^'s,  we  have 


I 


1a3xl 

Zap^lcsule 


L(s«oj  —  s)  C2  ^  a^e  C^3 
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(4.31) 


4  *  Za^ajcwe  t(als  “  C2  “  9  V  * 


Following  the  substitution  of  TTf  (4.29)  with  the  aid  of  Equa¬ 
tions  (4.23),  (4.28),  and  (4.31)  into  Equations  (3.45)  through  (3.50), 
we  finally  obtain  the  thermal  stress  field  as  follows. 
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*23  _  “wX3  -SX,i 
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*31  —  "o»X^  ”SX-1 
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It  can  be  shown  that,  for  a  single  moving  asperity,  i.e.,  N  *  1, 


Equations  (4.32)  check  with  the  Equations  (3.67)  for  the  case  =  0. 


4.5  Numerical  Results 


As  in  Chapter  3.0  the  mechanical,  the  temperature,  and  the  thermal 
solutions  are  all  expressed  in  the  Fourier  transformed  expressions. 
Inverse  Fourier  transform  will  be  applied  accordingly.  Numerical 
evaluation  of  the  integrals  and  corresponding  graphical  results  for 
Cases  I,  2,  and  3  will  be  presented. 


4.5.1  Rectangular  Contact  Area  with  Uniform  Pressure  (Case  1) 
Mechanical  Stress  Field 

In  view  of  Equation  (4.1),  the  transformed  boundary  conditions  are 


( j-l)+a 
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dx 


1 


2P. 


*  _  sinlx^a)  sin(x2b) 


WXjX2 


N  lx.x  (j-D 
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Rl  ’  "uf  R3  * 


(4.33) 


Substituting  the  expressions  above  into  Equations  (4.2),  applying  the 
Inverse  transform,  and  going  through  the  same  simplification  as  shown 
in  Section  (3.5.1),  leads  to  the  following  expression. 
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where 
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{[a  cos  $  +  bsin$  +  Y  cos(8  -  p)  -  x  (j 

Et  ■  In - 1 — 
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- - - ytr- 
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*,«“1  a  cos  8  +  b  sin  $  +  y  cos(8  -  6)  +  \  (j  -  1)  cos 


Temperature  Field.  The  expression  of  In  Equation  (4.20)  Is  equlva 

lent  to  -UfV  Is.  i.e., 

“2wfPnV  N  ix.x' (j-1) 

3  — 33 —  s1n(x1a)  s1n(x2b)  •  £  e  •  (4.3 

^x^2  j3l 


Substituting  Equation  (4.35)  into  (4.20)  and  taking  the  inverse  trans 
form  leads  to 
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where 
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M2  a  - 2  8  cos  x1Cx1  -  x  (j  -  1)]  +  sin  x1[x1  -  x  (j  -  1)]  , 

By  analogy  with  Equation  (3.72),  the  temperature  field  from 
Green's  function  approach  is 
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Thermal  Stress  Field.  After  substituting  the  expression  TJn  into 
Equation  (4.32)  and  applying  the  Inverse  transform,  we  find 
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4.5.2  Disk  Contact  Area  with  Uniform  Pressure  (Case  2) 
Mechanical  Stress  Field.  In  this  case  the  transformed  boundary 
conditions  are 
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(4.39) 


With  the  expressions  (4.39)  in  (4.2)  and  by  taking  Inverse  transform, 
the  solution  for  the  mechanical  stress  field  becomes 
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Temperature  Field.  With  the  expression  of  "JJn  given  by 
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(4.41) 


the  temperature  solution  can  be  readily  obtained  as 
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In  view  of  Equation  (4.12),  the  Green's  function  solution  Is 
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Thermal  Stress  Field.  The  expressions  are 
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4.5.3  Disk  Contact  Area  with  Paraboloidal  Distribution  of  Pressure 

VTsTS) - — 

Mechanical  Stress  Field.  Similarly,  we  can  obtain  the  transformed 
expressions  of  boundary  conditions  In  Equation  (4.1).  This  Is  given  by 
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(4.45) 


The  mechanical  stresses  are  analogous  to  those  outlined  in  Equa- 

ii  u  a 

tlons  (4.40)  except  and  It  are  replaced  by  Et  and  It,  respec¬ 
tively.  Their  forms  are 
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Temperature  Field.  According  to  Equation  (4.45),  we  have 
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Thus  the  temperature  field  with  the  aid  of  the  above  equation  after  the 


inverse  transform  is 


-2W4:P  V  N 
Hf  o  v 


ttIT 


N  /•»  /•< 

L  f  f 

3-iJo  Jo 


f3  cos(72x2)  [M^  cos(qx3) 


+  M2  sin(ex3)]  dx1  dx2  . 


(4.47) 


The  Green's  function  solution  Is 
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Thermal  Stress  Field.  Replacing  f2  in  Equation  (4.44)  with  f3  leads  to 
the  solution  accordingly. 


As  above,  the  temperature  and  the  thermal  stress  fields  for  three 
cases  have  the  form  of  the  summation  for  the  index  j  from  1  to  N  before 
double  Integrals,  which  would  Increase  computing  time  and  cost  for 
large  N.  This  Inefficiency  can  be  overcome  by  using  the  fact  that 
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where  tn  *  1/2  x' (N  -  1)  7X.  Thus  the  expressions  (4.35),  (4.41), 
and  (4.46)  for  ^  can  be  modified.  As  a  result,  the  temperature  and 
the  thermal  stress  fields  can  be  rewritten  as 
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where  r  *  1,  2,  and  3  for  Cases  1,  2,  and  3,  respectively.  Ai 
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For  the  numerical  results  which  follow  It  Is  considered  that  con¬ 
tact  Is  In  the  form  of  a  sequence  of  three  asperities  with  centers 
separated  by  a  distance  \  3  12.  The  values  used  for  the  moving 
velocity  of  asperities,  the  dimension  of  contact  areas,  and  the  coeffi¬ 
cient  of  friction  are  the  same  as  in  previous  chapters,  i.e.,  V  3  15 
m/s  {600  ips),  a  3  0.25  mm  (0.01  in.),  and  3  0.5,  respectively. 
Figure  4.2  shows  the  mechanical  principal  stress  fields  for  three  cases 
at  the  depth  c  3  10”1,  which  is  of  interest  here.  It  can  be  seen 
that  the  tensile  stresses  occur  on  the  trailing  side  and  the  compres¬ 
sive  stresses  on  the  leading  side  of  each  moving  asperity.  In  each 
case  the  corresponding  value  at  the  corresponding  location  beneath  each 
asperity  has  the  same  order  of  magnitude,  which  justifies  that  the 
cumulative  mechanical  effect  is  negligible  as  expected.  In  Figure  4.3 
is  shown  the  temperature  field  for  Case  1  in  the  surface  layer.  It  is 
noted  that  the  temperature  is  high  near  the  surface,  however,  its 
magnitude  drops  rapidly  with  respect  to  the  depth.  The  temperature 
gradient,  which  affects  the  thermal  stress  field,  attains  its  maximum 
at  c  3  10“ l.  Consequently,  high  tensile  stresses  are  again  expected 
at  this  depth.  The  temperature  rise  underneath  the  second  and  third 
asperities  is  due  to  the  cumulative  effect  from  the  previous  asperity. 
For  comparison,  the  temperature  fields  for  three  cases  are  presented  in 
Figure  4.4.  It  shows  that  Case  3  for  the  non-uniform  heating  has 
higher  temperature  than  Cases  1  and  2  for  the  uniform  heating.  In  all 
cases,  the  temperature  rises  at  the  corresponding  location  for  each 
asperity.  As  a  result,  large  tensile  stresses  occur  near  the  first 
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asperity  and  become  smaller  when  5  Is  large  (see  Figure  4.5),  The 
combined  effects  from  mechanical  and  thermal  stress  fields  for  three 
cases  are  shown  in  Figure  4.6.  Accordingly,  the  maximum  principal 
stress  for  Cases  1,  2,  and  3  would  reach  the  ultimate  strength  at  cri¬ 
tical  pressures  of  119  MPa  (17,245  psi),  104  MPa  (15,091  psl),  and 
88  MPa  (12,700  psl),  respectively. 


1  stress 


I 


[I 


LI 


5.1  Asperity-Size,  Shape  of  Contact  Area,  Load  Distribution 

Of  interest  here  are  the  effects  of  different  sizes,  shapes  of 
contact  areas,  and  load  distributions  of  a  moving  asperity  in  addition 
to  what  we  have  discussed  in  Section  2.5  (see  Figure  2.8  and  2.9).  In 
three-dimensional  analysis,  we  recall  that  the  pressure  distribution  Is 
non-uniform  for  Case  3  and  uniform  for  Cases  1  and  2;  while  the  shape 
of  the  contact  area  is  rectangular  for  Case  1  and  circular  for  Cases  2 
and  3.  For  a  comparison  among  the  three  cases.  Figures  5.1  and  5.3 
show  maximum  mechanical  and  thermal  principal  stress  fields,  respec¬ 
tively.  Moreover,  temperature  fields  are  presented  in  Figure  5.2  and 
maximum  combined  thermomechanical  principal  stress  fields  In  Figure 
5.4.  Accordingly,  we  observe  that  the  order  of  magnitude  of  the  solu¬ 
tions  from  Case  2  Is  close  to  (a  little  higher  than)  that  from  Case  1; 
while  the  solutions  from  Case  3  have  much  higher  order  of  magnitude 
compared  with  Cases  1  and  2.  Consequently,  it  can  be  argued  that  the 
shape  of  contact  area  Is  less  of  a  deciding  factor  than  the  pressure 
distribution  over  the  contact  area.  Since  the  inevitable  Irregulari¬ 
ties  in  the  contact  surfaces  will  cause  the  load  distribution  to  be 
non-uniform  when  two  sliding  solids  are  nominally  in  contact  over  some 
small  areas,  this  result  Is  of  great  significance  to  account  for  the 
strong  tendency  of  the  crack  onset  in  actual  practice.  In  order  to 
show  the  asperity  size  effect,  temperature  and  maximum  combined  thermo- 
mechanical  principal  stress  fields  from  three  different  sizes  are 


sizes  are  Investigated  in  Figures  5.5  and  5.6,  respectively,  where 
*i's,  i  *  1,  2,  3,  are  the  half  widths  of  the  contact  area,  which  are 
0.25  nm  (0.01  in.),  0.5  im  (0.02  in.),  and  0.125  mm  (0.005  in.), 
respectively.  Mention  should  be  made  that  the  nondimensionalization 
for  (xltx2,x3)  and  T  are  slightly  different.  Here  we  define  (^.m.c^) 
as  (x1,x2,x3)/^  and  ^  as  T^k/q  a.j.  The  comparison  should  therefore 
be  made  cautiously.  Referring  to  Figure  5.5  where  the  constant  average 
pressure  P0  and  heat  flux  q0  are  used,  it  can  be  shown  that  at  the 
same  location  higher  temperature  results  from  larger  asperity  size  if 
we  compare  (Sx,^)  with  (51/2,2<j»2)  and  (2c x, <^3/2)  and  take  into  account 
the  fact  that  the  temperature  decreases  as  the  depth  increases.  In 
fact,  temperature  is  by  all  means  expected  to  be  higher  when  the  total 
heat  input  over  the  contact  region  increases.  It  is  clear  that  higher 
surface  temperature  causes  higher  compressive  thermal  stresses  in  the 
thermal  layer  and  as  a  result,  the  tensile  stresses  at  ^  *  10" 1  from 
larger  asperity  size  are  higher  as  shown  in  Figure  5.6.  In  other 
words,  the  thermal  fracturing  is  more  likely  to  happen  if  the  total 
heat  input  increase  is  achieved  by  the  increase  of  the  contact  area 
over  which  the  average  pressure  is  assumed  to  be  constant.  This  result 
checks  with  Figure  2.9  obtained  from  the  two-dimensional  analysis. 

5.2  Comparison  of  Two-  and  Three-Dimensional  Analyses 

In  the  two-dimensional  analysis,  the  aspect  ratio,  t  a  b/a,  of  the 
moving  asperity  has  a  length  much  larger  than  the  width  measured  in  the 
direction  of  motion,  that  is,  a  moving  line  asperity  is  considered. 
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Maximum  combined  thermomechanical  principal  stress  fields  from 
different  asperity  sizes. 


For  the  mechanical  effect  portion,  the  dynamic  stress  potentials  are  of 
the  Lamd  type.  When  transformed  to  the  coordinates  fixed  to  the  asper¬ 
ity,  the  stress  potentials  become  harmonic  with  distorted  depth  coor¬ 
dinates.  The  distortion  is  caused  by  the  subsonic  velocity  of  the  mov¬ 
ing  asperity.  Using  a  unit  normal  and  a  unit  tangent,  the  Impulse 
responses  in  the  stress  and  the  displacement  fields  were  determined  to 
be  of  the  Eringen-Suhubi  type  [19].  The  impulse  responses  were 
employed  as  Green's  function  to  compute  the  displacement  and  the  stress 
fields  for  any  distribution  of  pressure  and  friction  force.  In  the 
thermal  effect,  following  Ling's  assumption  [48]  for  high  Peclet  num¬ 
ber,  the  heat  conduction  in  the  traverse  direction  is  neglected, 
retaining  only  the  convective  term  in  that  direction  from  the  local 
variation.  The  temperature  and  the  thermal  stress  fields  are  computed 
from  the  thermal  Impulse  responses,  which  were  obtained  from  a  unit 
heat  input  [22]. 

In  the  three-dimensional  model  of  a  single  moving  asperity  the 
potential  theory  cannot  be  readily  used,  that  Is  the  three-dimensional 
formulation  cannot  be  simply  extended  from  the  two-dimensional  mathema¬ 
tical  model.  Double  Fourier  transform  methods  are  therefore  applied 
for  deriving  the  analytical  solutions.  Due  to  the  complexity  of  the 
formulations,  they  are  first  left  in  the  transformed  space.  Numerical 
integrations  are  then  to  be  carried  out  for  corresponding  inverse 
transforms. 

It  is  clear  that  two-dimensional  theory  provides  a  much  simpler 
approach  which  also  leads  to  fairly  good  numerical  results  of  the 
mechanical  and  the  temperature  fields.  This  can  be  seen  if  we  compare 


Figures  2.5  and  2.6  with  3.2  and  3.5,  respectively.  However,  even  with 
the  analytical  difficulty,  the  three-dimensional  modeling  is  virtually 
unavoidable  as  it  best  describes  the  geometrical  shape  of  the  asper¬ 
ity.  Therefore,  in  that  model,  different  load  distributions  over  a 
circular  or  a  rectangular  area,  instead  of  the  unrealistic  line  asper¬ 
ity,  can  be  analyzed.  In  addition,  taking  Case  1,  t  *  10  from  the 
three-dimensional  theory  and  the  one  with  uniform  pressure  distribution 
from  two-dimensional  theory  for  comparison,  thermal  stresses  have  a 
different  order  of  magnitude  as  shown  in  Figures  3.8  and  2.8.  This  is 
deduced  from  the  dimensional  effect  which  is  of  great  importance.  As 
we  know  there  is  no  temperature  variation  in  the  n  direction  in  two- 
dimensional  analysis;  while  in  three-dimensions  the  temperature  drops 
off  near  the  edge  of  the  contact  area  in  the  n  direction,  i.e.,  n  * 

±10,  illustrated  by  Figures  5.7  and  5.8  for  Case  1,  t  *  10  at  depths  z 
«  10" 6  and  10" l,  respectively.  Thus  far,  the  cool  regions  outside 
the  hot  spot  in  the  n  direction  also  constrain  the  thermal  layer  from 
expanding.  As  a  result, . tensile  stresses  at  5  =  10"1  is  elevated  to 
a  higher  order  of  magnitude.  Thus  the  fracture  threshold  predicted 
from  three-dimensional  theory  is  at  a  much  lower  asperity  pressure  or  a 
lower  asperity  rubbing  speed  than  that  from  the  two-dimensional  theory. 
It  is  worth  mentioning  that  three-dimensional  theory  can  be  applied  to 
justify  two-dimensional  theory.  Numerically  it  is  shown  in  Figure  5.9 
where  the  result  for  Case  1  "moves  toward"  two-dimensional  solution  as 
the  aspect  ratio  increases.  Theoretically  the  Dirac  function  6(72) 
will  be  associated  with  the  transformed  boundary  functions  such  as  7i 
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in  Equation  (3.15)  or  Pi  in  Equation  (3.67)  if  we  consider  a  moving 
line  asperity  instead  of  a  moving  asperity  with  a  rectangular  contact 
area.  Consequently,  three-dimensional  formulations  of  double  integrals 
can  be  easily  simplified  to  single  integrals  which  converge  to  the  two- 
dimensional  solutions. 

5.3  Cumulative  Effects  of  Periodic  Loads 

It  is  evident  that  the  spacing  between  asperities  and  the  number 
of  repetitions  serve  as  governing  parameters.  The  cumulative  effects 
on  mechanical  stresses,  temperature  distributions,  and  thermal  stresses 
have  been  demonstrated  with  a  sequence  of  three  asperities  with  center 
distance  x'  3  12  in  Figures  4.2,  4.4,  and  4.5,  respectively.  In 
summary,  the  mechanical  stress  distribution  duplicates  individual 
effect  as  in  the  problem  of  a  single  moving  asperity.  The  temperature 
field  however,  because  of  the  periodic  updating  of  the  initial  tempera¬ 
ture  field  at  individual  asperities,  shows  a  definite  cumulative 
effect,  albeit  with  a  reducing  increment  at  individual  passing 
asperities.  Since  thermal  stress  distribution  is  sensitive  to  the 
temperature  variation,  the  magnitudes  of  tensile  stresses  in  the  vicin¬ 
ity  of  the  first  moving  asperity  are  increased.  As  a  result,  the  com¬ 
bined  stress  field,  as  shown  in  Figure  4.6,  reaches  a  fracture  thres¬ 
hold  at  a  much  lower  critical  average  pressure  compared  with  that  for  a 
single  moving  asperity.  If  the  center  distance  of  asperities  is 
increased  to  1200  times  the  characteristic  dimension  of  the  contact 
area,  i.e.,  x'  3  1200,  not  only  the  mechanical  stress  field  from 


individual  excitation  propagating  at  the  sound  speed  is  noncumulatlve, 
but  also  the  thermal  responses  from  the  preceding  asperities  have  died 
out.  Hence  In  this  case,  the  mechanical  and  thermal  responses  are 
essentially  repeated  under  energy  passage  of  moving  asperities.  As  for 
the  cumulative  effects  from  the  number  of  repetitions.  Figures  5.10  and 
5.11  show  temperature  fields  and  maximisn  combined  thermomechanical 
principal  stress  fields  for  different  number  of  moving  asperities.  As 
expected,  the  increase  of  the  number  of  repetitions  may  result  in 
higher  tensile  stresses  and  thus  the  seal  is  most  likely  easier  to 
crack.  Nevertheless,  the  damaging  aspect  of  cumulative  effects  Is  not 
to  Increase  Indefinitely.  For  If  we  recall  the  temperature  field  in 
Equation  (4.20) 

Qu  "uX* 

T  »  e  [C1  cos(ex3)  +  C2  s1n(ex3)] 
where  On  Is  associated  with  a  series  of  the  form 

N  ix.X  (j-1) 

E  e  1 

j-1 

which  has  a  close  form  as 

e  N  •  sindj-  NX  x^j/sln^  x  x^) 

with  tn  expressed  as  1/2  x’(N  -  1)  71#  It  can  be  shown  that  all 
terms  in  7  are  bounded  and  the  resultant  thermal  solutions  after  the 
inverse  Fourier  transform  are  convergent  on  the  consideration  of  the 
Riemann-Lebesgue  theorem,  as  the  number  of  periodic  excitation  of 
traversing  asperities  approaches  infinity. 
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Maximum  combined  therroomechanl cal  principal  stress  fields  due  to 
different  number  of  moving  asperities. 


Thus  far,  three  models  characterizing  the  heat  checking  failure 


have  been  established.  These  models  contribute  general  straight¬ 
forward  methods  so  that  we  can  assess  the  tendencies  of  the  crack  ini¬ 
tiation,  which  Is  of  major  concern.  The  avoidance  of  thermocracks  by 
controlling  some  parameters  such  as  E,  a,  k,  V,  a,  etc.  can  be  opti¬ 
mized.  In  general.  It  can  be  concluded  that  reducing  the  size  and  num¬ 
ber  of  asperities,  decreasing  the  moving  velocity,  and  having  uniform 
pressure  distributions,  and  choosing  a  seal  material  with  low  modulus 
of  elasticity,  low  coefficient  of  thermal  expansion,  high  thermal  con¬ 
ductivity,  and  high  fracture  strength  are  efficient  to  alleviate  the 
thermal  fracturing  problem. 

5.4  Justifications  and  Recommendations 

Of  concern  here  are  the  justifications  of  the  major  assumptions, 
the  formulations,  and  the  numerical  results.  In  the  analysis,  the 
flash  temperature  in  the  neighborhood  of  the  asperities  contained  on 
the  mating  ring  is  so  high  that  all  the  frictional  heat  may  be  assumed 
to  flow  predominately  into  the  primary  ring  since  the  temperature  dif¬ 
ference  between  the  asperities  and  the  primary  ring  is  much  higher  than 
that  between  the  asperities  and  the  mating  ring.  As  a  result,  this 
yields  an  upper  limit  of  the  solution  and  gives  a  conservative  estimate 
of  the  critical  pressure  which  would  result  in  heat  checking. 

In  practice,  the  velocity  of  the  moving  asperities  is  much  less 
than  that  of  the  Rayleigh  waves,  and  there  is  no  plastic  deformation 
observed  in  tension  tests  at  temperatures  as  high  as  1960#R  (1089  K) 
for  Stellite  III.  Consequently,  the  physical  conditions  of  the  low 


dilatation  rate  and  no  thermal  shock  phenomenon  justify  the  application 
of  the  theory  of  uncoupled  quasi-static  thermoelasticity.  In  addition, 
the  linear  theory  holds  since  the  base  material  subjected  to  the  asper¬ 
ity  friction  Is  essentially  elastic  from  observations  of  failed  speci¬ 
mens. 

In  order  to  justify  the  analytical  and  the  numerical  results 
obtained,  several  strategies  were  used.  Both  the  mechanical  and  the 
thermal  solutions  for  Case  1,  with  the  aspect  ratio  t  •*>  from  the 
three-dimensional  theory  were  checked  by  the  solutions  from  the  two- 
dimensional  theory.  The  general  solution  and  the  Integral  solution  of 
the  mechanical  stress  field  for  Case  1  from  a  single  moving  asperity 
were  checked  by  Eason's  result  [21].  He  fixed  the  coordinate  system  to 
the  semi -Infinite  solid  and  used  the  technique  of  the  triple  Fourier 
transform  to  analyze  the  mechanical  responses  from  a  moving  surface 
force.  As  for  the  temperature  field,  using  the  Green's  function  con¬ 
cept  or  the  Fourier  transform  technique  led  to  the  same  result.  In  the 
computation  of  the  thermal  stress  field  of  a  single  moving  asperity, 
the  approach  adopted  a  semi-inverse  scheme  by  assuming  the  displacement 
expression  to  satisfy  the  governing  equations,  the  boundary  and  the 
regularity  conditions.  The  results  were  checked  by  the  solutions  from 
periodic  excitation  with  N*1  using  the  method  of  the  Inverse  Differen¬ 
tial  Operator  in  conjunction  with  the  Cramer's  rule. 

The  numerical  results  were  also  carefully  obtained  and  justified. 
In  general,  the  numerical  evaluations  were  checked  with  both  Simp  and 
Quanc8  codes  for  single  integrals  and  with  both  the  double  Iterative 


Gausslan-Laguerre  quadratures  and  the  double  adaptive  Simpson  schemes 
for  double  Integrals  (see  Appendix  II). 

The  proposed  models  have  been  successfully  applied  to  study  the 
brittle  failure  In  the  submarine  face  seals.  The  material  properties 
are  considered  to  be  constant  with  respect  to  the  temperature  in  the 
analysis.  This  holds  since  the  temperature  above  the  stress-free 
ambient  temperature  at  the  depth  of  the  order  of  one- twentieth  asperity 
size  where  the  cracking  occurs  is  much  lower  than  that  In  the  thermal 
layer.  However,  it  may  be  of  practical  interest  to  Investigate  the 
multilayer  effect.  In  this  case,  the  material  properties  may  be 
constant  in  each  layer  but  varied  In  different  layers.  Mention  should 
also  be  made  that  the  analyses  so  far  are  restricted  to  the  homogeneous 
material.  Since  It  Is  very  hard  to  perfect  the  material  beneath  the 
surface,  it  is  worthwhile  to  study  the  layered  media  with  the  hetergen- 
eity  induced  by  the  rigid  inclusion  or  the  cavity  among  layers.  The 
principal  stress  field  from  combined  thermomechanical  effects  strongly 
depends  on  the  pressure  distributions  over  the  contact  areas;  while  the 
pressure  distributions  are  much  affected  by  the  contact  force  and  the 
geometries  of  the  asperities.  For  the  mechanical  face  seal,  the  total 
force  remains  constant  during  the  operation.  Therefore,  the  experimen¬ 
tal  research  to  optimize  the  control  of  the  size  of  the  moving  asperi¬ 
ties  Is  extremely  important. 
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METHODS  OF  INTEGRAL  TRANSrQRMS  AND  GREEN'S  FUNCTIONS 


If  we  have  a  function 
transform  Is  defined  by 


the  n-dimenslonal  Fourier 


expCi (x  •  x)]  dx 


and  the  appropriate  Inversion  form  Is 


where 


f(x)  * — i — -  f  f  7(7)  exp[-i(x  .  x)]  dx 
~  {/t T)"-'— ~  ~ 

*" 


(£.£)■  X1X1  <•  XjX2  *■  ...  ♦  X„X„ 


denotes  the  Inner  product  of  the  vectors  £  and  x  of 

In  applications  to  physical  problems  the  following  two  formulas 
are  often  employed  with  the  transform. 


D*  f(x)  »  (-I)M  7°  7(7)  , 

x  ~  ~  ~ 


xa  f(x)  *  (-1 ) Ia(  Da  7(7) 

7 

where  the  notations  are  explained  as  follows.  For  nonnegative  integers 
«jU  •  1.2 . n), 

a  a  ( t  •  •  • » an ) 

is  called  a  multi -Index  of  n  variables,  and 


|o|  3  aj  +  <*2  +  •••  +  an 

is  called  its  length. 

When  x  3  {x1,x2,...»xn)  is  a  point  in  ^n,  let  ;c°  be  the  monomial 
of  degree  |aj,  which  is 


X  -  xL  x2  ...  xn 

and  D“  be  the  differential  operator  of  order  jaj  with  respect  to  x. 


which  is 


For  example,  if  m  is  a  natural  number,  and  aa  are  constants  with 
indices  a,  then 

E  «  Da 

|a|<m  a 

expresses  a  partial  differential  operator  of  order  m  with  constant 
coefficients. 

Detailed  theories  and  applications  of  Fourier  Integral  transform 
have  been  well  established  by  Snedden  [49]  and  many  other  authors.  The 
aforementioned  formulation  serve  as  the  basic  tool  and  are  sufficient 
for  the  analysis. 

The  Green's  function  used  to  solve  the  partial  differential  equa¬ 
tions  of  mathematical  physics  can  be  found  in  several  references  [50- 
53],  Of  concern  here  is  the  use  of  Green's  function  In  the  solution  of 
the  general  heat  conduction  equation.  For  generality  the  problem  is 
formulated  as 
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7  T(r,t)  +  ^  g(r,t)  *  -  — ^ —  In  R»  t  >  0  , 

k1  %+  hiT  “  f1(I*t}  on  S1*  t  >  0  » 

T(r,t)  *  F(r)  for  t  *  0,  In  R 

where  3/an-j  denotes  differentiation  along  the  outward-drawn  normal  to 
the  boundary  surface  S,,  i  *  l,2,..-.,m  and  m  is  the  number  of  contin¬ 
uous  boundary  surfaces  of  the  region  R.  g(jr,t),  the  internal  heat 
generation,  and  ff(£,t),  the  boundary-condition  functon  on  S-j,  vary 
with  position  and  time.  F <jr )  is  the  initial  condition  function. 

The  corresponding  Green's  function  is  the  solution  of  the  follow¬ 
ing  auxiliary  problem  for  the  same  region  R: 

+  i  6(r  -  r')  «(t  -  t)  *  In  R,  t  >  T  , 

ki  l§7 *  V5  ■ 0  *" 

G  *  0  If  t  <  t 

where 

I  III 

<5(r  “  X  )  *  SUj  -  XjJ  S(x2  -  x2)  $(x3  -  x3) 
for  (Xi.x2.x3)  coordinates. 


It  Is  noted  that  this  auxiliary  problem  satisfied  by  the  Green's 
function  has  homogeneous  boundary  conditions  and  has  an  impulsive  point 
heat  source  and  zero  initial  condition.  In  other  words,  this  solution 


represents  the  temperature  distribution  in  the  region  R,  which  is 
initially  at  zero  temperature  and  subjected  to  homogeneous  boundary 
conditions,  due  to  an  impulsive  heat  source  of  strength  unity  located 
at  jr*  and  releasing  its  heat  spontaneously  at  time  t. 

The  temperature  distribution  in  terms  of  the  Green's  function 
according  to  Ozisik  [54],  is 

T(r,t)  *  f  G|  _n  F(r  )  dv  +  f  dr  f  Gg(r  ,t)  dv 
IT-0  ~  k  /  n  J „  3  ~ 


rtm 

/  dt  E  f  fi,  f  Mr.t)  dS. 
J  T*0  1  =  1  Js<  Si  *1  1  ~  1 


where  G j refers  to  the  Green's  function  evaluated  at  the  boundary 
surface. 

Most  often  the  heat  conduction  equation  is  not  in  the  form  pre¬ 
sented.  Under  this  circumstance,  a  different  Green's  function  needs  to 
be  determined  or  some  transformation  may  be  useful  to  reduce  the  equa¬ 
tion  to  the  standard  form  for  which  Green's  function  is  known.  For 
example,  if  we  have 

at  *7  1  si  axx  h  ax2  03  ax3  Y'  9 

where  k,  8p  82*  83  and  T  are  constants.  We  define  a  new  dependent 
variable  u  as 


'm*  •**•*/*  m*m  •'*  .*•  . 


APPENDIX  II 


NUMERICAL  INTEGRATION 


The  numerical  methods  for  evaluating  the  integral  of  the  form 


I 


f (x)  dx 


9 


with  [a,b]  finite,  are  well  established  and  can  be  found  in  most  text¬ 
books  of  numerical  analysis  [55-59],  In  general,  the  Integration  meth¬ 
ods  can  be  viewed  within  a  simple  framework.  First,  try  to  find  a 
family  of  functions,  {fn(x),n  j>  1},  of  which  each  approximates  the 
Integrand  f(x)  in  the  following  sense 

if  -  ♦  0  as  n  ♦  •  . 


The  form  of  each  fn  should  be  chosen  such  that  we  can  easily  evaluate 
the  Integral 

In(f)  * J  fn(x)  dx 
a 

which  approximates  1(f),  that  Is  1(f)  •  In(f).  As  for  the  error,  we 
have 


En(f)  =  1(f)  -  In(f)  *  J  [f ( x)  -  fn(x)]  dx 

cl 


-  fn(x)  I  dx  <  (b  -  a)  if  -  f n . 


and 


The  Epical  best-known  example  is  the  Newton-Cotes  Integration 
formulas  which  are  obtained  by  a  suitable  interpolation  polynomial 
fn(x)  on  the  nodes  x0,x1,...,xn  to  replace  the  integrand  f(x). 
Consider  an  equally-spaced  partition  of  the  closed  interval  [a,b] 
given  by 

xi  a  a  +  *  fi  *  *  s  0,1*. ..,n  . 


Let  fn  be  expressed  in  connection  with  the  Lagrange  polynomial  as 


f  n  (x)  *  t  f  (Xj  )  Lf(x)  ,  L.U)  a  77  . 

i*0  1  j*0  xi  xj 

j*i 


If  we  define 

In(f)  -  £  f(xi ) 

i»0  1  1 


and  introduce  the  new  variable  t  such  that  x  *  a  +  t(b  -  a),  then  a ^ 
can  be  expressed  as 

«4  ■  (b  -  a)  f1  n  dt 

•'0  j*0  J 

which  depends  on  n  only.  For  sufficiently  smooth  f(x),  it  can  be  shown 
[60]  that  the  approximation  error  has  the  form 

h  a  P  +1  (PJ 

En(f)  *  VH-2-)  f  <«>  •  ?  CIO) 


where  Pn  and  Kn  depend  only  on  n  but  not  on  the  Integrand  f. 

In  particular,  for  n  3  2,  the  Newton-Cotes  formula  becomes  the 
Simpson  rule  which  is 
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J  f (x)  dx  -  s-g-s.  [f (a)  +  4f(^i)  +  f(b)]  -  ^QQ  fV‘t/  U> 


The  corresponding  composite  rule  is  obtained  as  follows 


r  b 

I  f(x)  dx 
•'a 


N-l  /**1+1 

£  J  f(x)  dx 
1=0  \ 


N-l  rl 

23  hi  /  f(x*  +  h.,t)  dt 

i*0  1  •'0  1  1 


N-l  h.  /  h.\  N 

*  £0  r  iflxi}  +  4f(xi +  r) +  f(xi+i)  ■  | 


tfi  fl4)  <c,) 

■  jT0  2555 


where  h^  *  x1+1  -  x, .  For  the  case  of  equally-spaced  intervals,  l.e., 


h<  -  *  h  , 


the  above  formula  becomes 


/-b  .  (  N-l  ,  N-l 

/  f (x)  dx  =  £  <f(a)  +  4  £  f(a  +  (k  +  i)  h)  +  2  £  (a  +  kh)  +  f(b) 

■'a  k»0  k=l 


(b  -  a)°  r(4) 


2880  N 


f'  '  (n)  ,  n  C(a.b)  . 


Another  powerful  Integration  method  Is  the  Gaussian  quadrature 
behind  which  the  Idea  is  to  find  an  Integration  formula 


In(f)  *  23  A.f(xJ  -  f  o»(x)  f(x)  dx 
j=l  J  J  •'a 


where  <u(x)  is  a  given  non-negative  weight  function  on  the  Interval 
(a,b),  which  may  be  Infinite;  and  assumes  the  following  properties: 


1.  I  1  x [ n  dj(x)  dx  is  Integrable  and  finite  for  all  n  ^  0. 
•'a 


2. 


for  some  continuous  non-negative 


function  g(x),  then  g(x)  =  0  on  (a,b). 

Aj,  Xj  can  be  determined  by  the  method  of  undetermined  coefficients 
which  asks  for  perfect  accuracy  when  f(x)  Is  one  of  the  power  functions 
l»x,x2,...,x2n”1.  This  provides  2n  conditions  for  the  determined  2n 
unknowns  Aj  and  Xj.  In  fact. 


aj 


Lj  (x)  dx 


where  Lj(x)  Is  the  Lagrange  polynomial.  The  arguments  x1,...,xn  are 
the  zeros  of  nth  degree  orthogonal  polynomials  *n(x)  on  (a,b)  with 
respect  to  the  weight  function,  l.e., 

ui(x)  a  (x)  *„,(x)  dx  »  0  for  n  #  m  . 

Tn  Tm 

The  error  expression  can  be  shown  to  have  the  form  [61] 


E  (f)  — f(2n)  (n)  ,  n  C(a.b)  . 

n  A‘(2n)l 


Gaussian-Laguerre  quadrature  is  one  particular  type  of  Gaussian  quadra 


ture  with  a»(x)  *  e"x.  It  takes  the  form 


the  arguments  x-j  being  the  zeros  of  the  nth  Laguerre  polynomial 


Pn(x)  *  eX  A-  <e"xxn) 
n  dxn 


and  the  coefficients  A*  being 


A 


1 


(nl)2 

xf[p;<xfn2  • 


The  numbers  xj  and  Ai  are  available  In  Tables  [62]. 

It  Is  known  that  an  Iterative  scheme  or  an  adaptive  scheme  should 
be  employed  In  conjunction  with  the  composite  rule  of  integration.  In 
an  iterative  scheme,  successive  approximations  to  the  Integral  are 
computed  until  It  has  met  the  required  accuracy.  For  Instance,  the 
algorithm  can  be  set  as  follows. 

Let  I2i<N(f)  denote  the  approximation  solution  based  on  the  compo¬ 
site  rule  with  2kN  equal-width  subintervals.  For  k  ■  0,1,2,...,  If 


1  k+1  (f)  ’  1  k  (f) 
2* 2kN 


is  sufficiently  small,  then  we  exit  with  I2k+lN^  as  the  result. 

The  Iterative  scheme  can  be  quite  inefficient,  since  the  accuracy 
may  be  satisfactory  on  many  subintervals  where  the  refinement  becomes 
unnecessary.  Based  on  this  Idea,  an  adaptive  scheme,  which  refines 
only  where  necessary,  has  been  well  established  [63-69].  The  Integra¬ 
tion  codes  Simp  and  Quanc8  associated  with  the  Simpson  rule  and  8  panel 
Newton-Cotes  formula,  respectively,  are  examples  of  such  a  scheme. 


Singular  Integrals,  whether  for  singular  Integrands  or  for 
infinite  range  of  Integration,  are  encountered  in  previous  chapters. 
Since  strategies  to  handle  such  singularities  are  too  varied  and 
numerous,  only  some  of  them  which  have  been  applied  successfully  are 
delineated  in  the  following. 

1.  Apparent  (or  removable)  singularities  may  exist  at  some  values 
of  the  integration  variable.  Limiting  values  should  be  found  with  the 
help  of  L'Hospital's  rule  if  necessary. 

2.  A  change  of  variable  may  exchange  a  difficult  singularity  for 
a  less  critical  one,  or  it  may  remove  the  singularity  completely.  For 
instance,  if  f(x)  £  C[0,1],  the  change  of  variable  tn  »  x  transforms 


the  integral 


.1  -I 


f  x  n  f (x)  dx  ,  n  >  2  , 
Jn 


n  f  f(tn)  tn"2  dt  , 
J  n 


which  Is  a  proper  integral.  Note  that  a  change  of  variable  Is  one  of 
the  most  effective  methods  to  evaluate  an  Integral  with  an  Infinite 
limit  since  it  can  transform  the  infinite  r  ion  into  a  finite  region. 


For  example,  if  y  ■  e-x,  then 


f(x)  dx 


f(-lnv) 


■'V/'-v’-.-'l-'- 


*.*  *_•  .*  V  V  r  / 


3.  Reduce  the  Infinite  Integral  fQm  f(x)  dx  to  a  finite  Integral 
f(x)  dx  and  estimate  the  truncation  error  due  to  f"  f(x)  dx  if 

possible. 

4.  It  has  been  suggested  by  Squire  [70]  that  the  Infinite  Inte¬ 
gral  can  be  transformed  to 

f  f (x)  dx  *  s  f  [ f ( st)  ♦  t"2  f ( s/t)]  dt  . 

*'0  •'o 


For  the  double  Integral  of  the  form 


f  f  f(x,y)  dx  dy  , 
•'0  •'0 


both  the  Inner  and  outer  Integrals  can  be  evaluated  by  using  Gausslan- 
Laguerre  quadrature  bas.ed  on  an  Iterative  scheme.  Alternatively,  the 
Integral  above  can  be  first  transformed  Into  the  Integral  with  finite 
Intervals  by  the  technique  of  change  of  variable.  Then  the  application 
of  any  of  the  adaptive  quadrature  schemes  of  two  dimensions  can  lead  to 
respectable  results. 
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